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Dynamic  aspects  of  a conceptual  Gas  Core  Reactor  (GCR)  are 

investigated  in  this  dissertation.  The  reactor  is  a bimodal,  closed 

cycle,  circulating  fuel  system  capable  of  operating  in  low  and  high 

power  modes.  It  uses  uranium  hexafluoride  (UFe)  in  gaseous  form  as 

the  fuel,  and  the  fuel  gas  mixture  of  UFe  and  helium  also  serves  as  the 
coolant /working  fluid. 

The  gaseous  nature  of  the  fuel,  plus  the  fact  that  this  is  a 
circulating  fuel  reactor,  makes  the  system  unique,  and  the  dynamic 
behavior  tends  to  be  quite  different  from  that  of  conventional,  solid 
core  systems.  The  system  behavior  while  undergoing  rapid  transition 
from  a low  power  mode  of  operation  to  a high  power  mode  of 
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operation  is  investigated  with  the  help  of  a computer  program 
developed  for  that  purpose.  The  program,  DYNAM,  written  in  standard 
FORTRAN  77,  simultaneously  solves  the  relevant  point  reactor 
kinetics,  thermod)mamic,  heat  transfer,  and  1-D  isentropic  flow 
equations  of  the  OCR. 

A nonlinear  stability  analysis  of  the  system  is  also  performed 
using  DYNAM.  The  results  of  this  analysis  indicate  that  the  investigated 
OCR  is  an  inherently  stable  system— mainly  due  to  the  strong  negative 
fuel  mass  feedback  effect.  Due  to  this  strong  fuel  mass  feedback,  the 
response  of  the  OCR  to  external  reactivity  perturbations  is  found  to  be 
quite  different  from  that  of  conventional,  solid  core  reactors.  For 
example,  at  full  power  operation,  it  is  observed  that  the  inherent 
reactor  response  is  a power  decrease  for  a positive,  as  well  as  a 
negative,  external  step  reactivity  insertion.  To  increase  GCR  power,  it 
is  necessary  to  provide  additional  positive  reactivity  during  the 
transient  by  varying  a system  operating  parameter  (for  example,  by 
increasing  the  core  inlet  pressure). 

A linear  stability  analysis  of  the  GCR  shows  that  the  linear 
analysis  is  often  adequate  for  correctly  predicting  system  stability. 
However,  for  certain  operating  regimes,  where  the  fuel  mass  feedback 
coefficient  is  large,  the  linear  analysis  breaks  down  due  to  the  strong 
nonlinear  effects.  For  these  regimes,  the  nonlinear  analysis  has  to  be 
performed  to  reliably  predict  the  system  stability. 
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CHAPTER  1 
INTRODUCTION 


Terrestrial  power  generation  making  use  of  nuclear  fission  is 
currently  accomplished  with  solid  fueled  reactors.  The  fuel  usually  is 
in  the  form  of  small  pellets  made  of  uranium  dioxide  (UO2)  or  in  some 

cases  uranium  metal  or  some  other  compounds  of  uranium  like 
uranium  carbide  or  nitride.  The  concept  of  using  nuclear  fuel  in 
gaseous  form  was  first  suggested  in  1955  (Ref.l).  Here  the  fuel  can  be 
either  uranium  plasma  (provided  the  temperature  in  the  system  is 
kept  high  enough  for  the  plasma  to  exist)  or  a compound  like  uranium 
hexafluoride  (UFe)  or  uranium  tetrafluoride  (UF4)  in  gaseous  form. 

Nuclear  fuel  in  gaseous  form  has  many  potential  advantages  over  that 
in  solid  form  which  could  eventually  make  it  the  fuel  form  of  choice 
for  many  special  applications  including  space  power. 

In  this  dissertation.  d)mamic  aspects  of  a conceptual  gas  core 
reactor  system  are  investigated.  The  reactor  uses  highly  enriched  UFe 

or  UF4  in  gaseous  form  as  the  fuel.  The  peak  fuel  temperatures 

attained  are  in  the  range  of  only  a few  thousands  of  degree  Kelvin  and, 
hence,  the  fuel  is  not  in  a plasma  state. 
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Knowledge  of  a reactor  system's  behavior  in  the  nonsteady  state 
is  of  primary  importance  in  the  design  and  construction  of  any  new 
type  of  reactor.  The  purpose  of  the  research  described  in  this 
dissertation  is  twofold.  The  first  goal  is  to  perform  the  system  startup 
analysis,  i.e.,  analysis  and  prediction  of  the  system  behavior  during 
various  startup  scenarios  and  the  second  goal  is  to  analyze  the  stability 
of  the  system's  high  power  mode  of  operation,  i.e.,  to  establish  the 
conditions  under  which  the  reactor  can  be  operated  in  a safe  and 
reliable  manner  at  high  power  levels. 

Advantages  of  Gaseous  Fuel  Over  Solid  Fuel 
One  of  the  major  limitations  of  conventional  solid  core 
pressurized  water  reactors  (PWR)  or  boiling  water  reactors  (BWR)  is 
their  relatively  low  working  fluid  temperature.  Since  the  zircaloy 
cladding  of  the  fuel  loses  its  strength  at  around  700  K,  the  working 
fluid/coolant  temperature  has  to  be  kept  below  this  value.  This  is  a low 
value  compared  to  the  1500  K or  so  reached  at  the  centerline  of  an 
average  fuel  pellet  inside  the  core,  and  this  results  in  a low  power 
plant  energy  conversion  efficiency  compared  to  a coal  or  oil  fired 
power  plant.  The  reason  for  not  replacing  zircaloy  with  some  other 
metal  or  alloy  with  better  mechanical  strength  is  that,  except  for  this 
drawback,  zircaloy  is  an  excellent  cladding  material  with  many 
desirable  properties  like  low  neutron  absorption  cross  section,  good 
resistance  to  corrosion  by  water  at  elevated  temperatures,  and 
reasonably  good  resistance  to  radiation  damage.  By  going  to  a gaseous 
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fuel,  this  serious  limitation  of  solid  core  systems  can  be  circumvented. 
A gaseous  core  system  can  operate  at  extremely  high  temperatures 
(and  high  efficiencies)  since  the  fuel,  either  by  itself  or  when  mixed 
with  another  gas,  can  also  be  the  coolant/ working  fluid,  and,  unlike  in 
a solid  core  system,  there  may  not  need  to  be  any  heat  transfer 
process  from  the  fuel  to  a coolant/working  fluid.  Practical  heat 
removal  and  materials  considerations  generally  force  one  to  operate 
the  gas  core  reactor  at  temperatures  in  the  range  of  a few  thousands  of 
degree  Kelvin.  An  exception  to  this  is  the  plasma  core  reactors  where 
the  fuel  temperatures  are  in  the  tens  of  thousands  of  degree  Kelvin 
range.  Here  special  arrangements  are  made  to  keep  the  fuel  away  from 
the  containment  wall. 

A gaseous  core  system  also  opens  up  the  possibility  of  continuous 
operation  of  the  power  plant  with  on-line  refuelling  and  on-line  fission 
product  removal.  Other  advantages  associated  with  gas  core  reactors 
include  a much  less  complicated  core  structure,  high  fuel  utilization, 
simple  fuel  management,  and  inherent  safety  due  to  the  negative 
reactivity  associated  with  an  expanding  gaseous  fuel. 

Special  Advantages  of  Gas  Core  Reactors  for 
Space  Power  Applications 

The  high  working  fluid  temperature  plus  the  gaseous  nature  of 
the  fuel  makes  the  system  quite  attractive  for  space  power 
applications.  First  of  all,  a high  working  fluid  temperature  implies  high 
efficiency,  which  in  turn  implies  a compact  system  for  a given  power 
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output.  A second  advantage  is  the  reduction  in  the  size  of  the  radiator. 
For  any  space  power  source,  nuclear  or  non-nuclear,  the  problem  of 
efficiently  rejecting  waste  heat  is  always  a major  concern.  Radiation  is 
the  only  way  of  rejecting  heat  in  space,  and,  for  large  power  systems, 
the  mass  of  the  radiator  usually  tends  to  dominate  the  total  system 
weight  (Ref. 2).  Thus,  there  is  a strong  incentive  for  reducing  the 
radiator  mass.  Since  the  rate  of  radiative  heat  transfer  varies  as  the 
fourth  power  of  temperature,  the  very  high  working  fluid 
temperatures  possible  with  gas  core  reactors  offer  significant 
advantages  with  regard  to  radiator  size  and,  hence,  system  mass. 

A third  advantage  is  the  elimination  of  the  accidental  criticality 
concerns.  Before  the  launching  of  any  reactor  into  space,  extreme 
precautions  have  to  be  taken  to  make  sure  that  the  system  remains 
subcritical  throughout  normal  launching  operations  as  well  as  during 
any  foreseeable  accidents  during  launch.  By  using  a gaseous  core 
system,  the  fuel  in  its  storage  container(s)  is  (are)  in  a highly 
subcritical  state  during  launch.  Achievement  of  criticality  requires  that 
the  gas  be  both  pressurized  and  contained  or  surrounded  by  an 
efficient  moderator/reflector.  Also,  since  the  process  of  putting  a 
reactor  in  orbit  is  an  extremely  costly  affair,  one  would  like  the  reactor 
to  remain  in  operation  for  a long  period  of  time.  The  gas  core  reactor 
offers  the  possibility  of  long  term  operation  without  building  in  a large 
amount  of  excess  reactivity  at  the  time  of  launch  due  to  its  capability 
for  on-line  refuelling  and  fission  product  removal. 


Present  Work  Description 

The  concept  analyzed  in  this  dissertation  is  a gaseous  core 
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reactor  power  system  that  is  designed  to  operate  at  two  widely 
different  power  levels:  a low  power  mode  of  operation  and  a burst  or 
high  power  mode  of  operation.  It  is  a highly  enriched  UFe  or  UF4 

fueled  reactor  capable  of  delivering  power  anywhere  from  a few  MWe 
up  to  a few  hundreds  of  MWe. 

The  advantage  of  using  UFe  is  that  it  is  the  only  stable  compound 

of  uranium  that  can  exist  in  a gaseous  state  at  a relatively  low 
temperature  (sublimation  point  = 329.5  K at  atmospheric  pressure). 
This  low  melting  point  implies  that  the  system  temperature,  including 
that  of  the  external  loop,  can  be  kept  at  a reasonably  low  value,  and 
still  avoid  condensation  of  the  fuel  on  the  inner  surfaces  of  the  pipes 
and  ducts.  The  limitations  of  UFe  are  that  it  is  highly  corrosive  and  it 

undergoes  significant  dissociation  by  2000  K.  As  for  UF4,  significant 

dissociation  starts  only  around  5000  K and  UF4  is  also  much  less 

corrosive  than  UFe . The  problem  in  using  UF4  as  the  fuel  is  that  it  is 

very  difficult  to  start  the  system  up  because  of  the  high  melting  and 
boiling  points  (1300  K and  1700  K,  respectively,  at  atmospheric 
pressure)  of  UF4.  One  possibility  is  that  the  system  could  be  started  up 

using  UFe  as  the  fuel  and  then  once  the  fuel  temperature  gets  high 

enough  (above  2000  K),  UF4  naturally  becomes  the  fuel  since  one  of 
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the  major  dissociation  products  of  UFe  is  UF4.  A more  detailed 

description  of  the  system  is  given  in  Chapter  2. 

The  incentive  for  designing  such  a bimodal  system  primarily 
came  from  Strategic  Defense  Initiative  (SDl)  requirements.  The  power 
requirements  specified  by  the  SDIO  can  be  divided  into  three  distinct 
ranges: 

a)  a station  keeping  mode  requiring  a few  MWe  for  a long  period 
of  time  (7-10  years). 

b)  an  alert  or  enhanced  surveillance  mode  requiring  a few  tens 
of  MWe  for  shorter  periods  of  time  (a  few  hours  to  a few  days)  and 
finally, 

c)  a burst  power  or  defense  mode  needing  a few  hundreds  of 
MWe  or  more  for  short  periods  of  time  (a  few  tens  of  minutes  up  to  an 
hour  or  so). 

The  bimodal  gas  core  reactor  system  that  is  described  in  this 
dissertation  is  expected  to  meet  all  of  the  above  power  requirements, 
i.e.,  the  low  power  mode  is  capable  of  meeting  the  power 
requirements  of  both  the  station  keeping  phase  as  well  as  the  alert 
phase  and  the  high  power  mode  of  operation  can  handle  the  power 
requirements  during  the  burst  power  phase. 

As  a part  of  the  system  design,  both  steady  state  (static)  analysis 
and  dynamic  analysis  have  to  be  performed.  This  dissertation  is  mainly 
concerned  with  the  dynamic  analysis  of  the  system. 
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Logically,  the  dynamic  analysis  can  be  divided  into  three  distinct 
phases: 

a)  analysis  of  the  system's  transient  behavior  during  the  low 
power  mode  of  operation. 

b)  analysis  of  the  transition  from  the  low  to  the  high  power 
mode  of  operation,  and  finally. 

c)  analysis  of  the  system's  transient  behavior  during  the  high 
power  mode  of  operation. 

Dynamic  analysis  of  the  low  power  mode  has  already  been 
performed  at  the  University  of  Florida  (Ref.3).  This  work,  therefore, 
concentrates  on  the  remaining  two  areas. 

Specifically,  the  first  part  of  the  dissertation  looks  at  the  startup 
of  the  high  power  (or  the  transition  to  the  high  power)  mode  of 
operation. 

One  of  the  requirements  which  any  space-based  power  system 
for  SDl  burst  power  applications  has  to  satisfy  is  the  capability  for 
rapid  startup,  i.e..  within  a short  period  of  time  the  system  should  be 
able  to  go  safely  and  in  a predictable  manner  from  a low  power, 
stationkeeping  mode  to  a high  power  mode.  SDIO  requirements  limit 
this  transition  time  to  about  100  seconds,  which  is  extremely  short 
compared  to  the  few  hours  taken  by  conventional,  terrestrial  nuclear 
plants  to  reach  full  power  operation  starting  from  low  power  or  hot 
standby  conditions. 
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The  response  of  the  gas  core  reactor  that  is  analyzed  in  this 
dissertation  is  expected  to  be  quite  different  from  that  of  a 
conventional,  solid  core  system.  The  highly  compressible  nature  of  the 
fuel,  along  with  the  fact  that  this  is  a circulating  fuel  system,  implies 
that  the  fuel  mass  within  the  core  is  not  going  to  be  a constant  during 
transients.  The  interaction  between  the  core  pressure,  the  core 
temperature,  and  the  fuel  mass  within  the  core  adds  additional 
non-linearity  to  the  system  equations.  It  is  very  important  to  know  the 
behavior  of  key  system  parameters  during  the  rapid  transition  from 
the  low  power  mode  to  the  high  power  mode  of  operation  in  order  to 
ensure  a safe  transition. 

To  analyze  the  dynamic  behavior  of  the  gaseous  core  nuclear 
reactor  during  rapid  transitions,  a computer  program  is  developed. 

The  program  (DYNAM)  solves  the  complete  system  equations 
(nonlinear)  numerically  and  predicts  the  behavior  of  important  system 
parameters  during  the  transition.  It  can  be  a valuable  tool  for  stud)nng 
the  effects  of  variations  in  different  system  parameters  on  the  system 
response  and  also  for  defining  safe  transition  paths. 

The  second  part  of  the  dissertation  is  concerned  with  the 
stability  analysis  of  the  high  power  mode  of  operation,  i.e.,  assuming 
that  a safe  transition  has  already  been  made  from  a low  power  mode  of 
operation  to  a high  power  mode  of  operation,  this  section  looks  at  the 
power  behavior  following  any  perturbations  (primarily,  reactivity 
perturbations)  imposed  on  the  steady  state,  high  power  system.  If  the 
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system  is  found  to  be  inherently  unstable,  the  system  design  has  to  be 
modified  to  make  it  stable. 

One  way  to  perform  the  stability  analysis  is  to  linearize  the 
nonlinear  system  equations,  and  then  use  the  powerful  methods  of 
linear  stability  analysis  to  determine  the  system  stability.  One  of  the 
objectives  of  this  section  is  to  determine  whether  a linear  stability 
analysis  is  adequate  or  not  for  predicting  the  stability  of  the  gaseous 
core  reactor  system.  To  this  end,  the  results  from  the  linear  analysis 
can  be  compared  against  those  obtained  using  DYNAM.  If  the  linear 
analysis  is  found  to  be  adequate,  one  can  use  it  to  determine  quickly 
the  range  of  values  for  reactor  physical  parameters  and  equilibrium 
power  levels  for  which  the  system  is  stable,  and  to  predict  the  bounds 
on  allowable  departures  from  the  equilibrium  states  without  having  to 
resort  to  nonlinear  analysis—which  can  get  very  involved  and  time 
consuming  for  a system  like  this. 

Dissertation  Organization 

Chapter  2 contains  a review  of  the  previous  work  done  on 
gaseous  core  reactor  power  systems.  It  includes  a review  of  the  gas 
core  reactor  work  that  has  been  going  on  at  the  University  of  Florida 
for  the  past  15  years  or  so.  It  also  includes  a detailed  description  of 
the  bimodal  reactor  that  is  the  topic  of  this  dissertation. 

This  is  followed  by  a chapter  which  describes  the  transition  or 
startup  analysis  of  the  high  power  system.  The  models  used  in  the 
computer  code  developed  for  analyzing  the  transition,  as  well  as  some 
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of  the  limitations  of  the  code,  are  discussed.  This  chapter  also 
presents  some  of  the  representative  results  generated  by  the  code. 

Chapter  4 contains  the  theory  and  a description  of  the  methods 
used  for  the  nonlinear  and  linear  stability  analysis  of  the  burst  power 
system.  It  also  contains  the  results  of  the  stability  analysis. 

The  last  chapter  presents  the  conclusions.  It  also  contains  some 
suggestions  for  future  work,  and  modifications  that  can  be  made  to  the 
computer  program,  DYNAM,  written  for  analyzing  the  gas  core  reactor. 


CHAPTER  2 

PREVIOUS  WORK  ON  GASEOUS  CORE  REACTORS 

Research  on  gaseous  core  reactors  has  included  many 
independent,  conceptual  phase  efforts  for  the  past  30  years  or  so. 
During  this  time,  a wide  variety  of  system  configurations  and  physical 
arrangements  for  the  reactor  core,  energy  conversion  schemes,  and 
fuel/working  fluid  form  have  been  suggested,  and  a wide  range  of 
scientific  feasibility  studies,  mostly  analytical  and  a few  experimental, 
have  been  conducted  with  encouraging  results.  As  mentioned  in 
Chapter  1,  in  the  United  States  the  gas  core  reactor  concept  was  first 
suggested  by  Bell  in  1955  (Ref.  1).  He  performed  preliminary  static 
neutronics  calculations  on  a spherical  core  fueled  with  UFe  gas  and 

surrounded  by  a moderator/reflector  made  of  either  heavy  water 
(D2O),  beryllium  or  graphite.  The  impetus  for  much  of  the  early  work 

on  gas  core  reactors  came  from  the  National  Aeronautics  and  Space 
Administration  (NASA).  Primarily,  NASA  was  interested  in  the 
possibility  of  applying  such  systems  for  rocket  propulsion  in  space. 

A number  of  concepts  were  evaluated  by  NASA  and  out  of  these, 
only  two  were  selected  for  further  study  (Ref.2).  These  two  were  the 
open-cycle  coaxial  flow  system  and  the  closed-cycle  nuclear  light  bulb 
concept.  The  general  idea  in  both  these  concepts  is  to  use  the  energy 
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from  a fissioning  plasma  (of  partially  ionized  uranium  with  average  fuel 
temperatures  in  excess  of  10,000  K)  to  radiatively  heat  the  propellant 
which  is  then  exhausted  through  a nozzle  to  generate  the  required 
thrust. 

The  original  work  on  the  open-cycle  coaxial  flow  reactor  was 
performed  by  Rom  (Ref. 4)  and  Ragsdale  (Ref. 5)  under  the  direction  of 
the  NASA  Lewis  Research  Center.  The  coaxial  flow  concept  involves  a 
spherical  cavity  with  a low  velocity  inner  stream  of  fissioning  uranium 
metal  plasma.  This  hot  plasma  is  kept  away  from  the  cavity  walls  by  a 
high  velocity  propellant  gas  (typically  hydrogen)  stream.  The  hot  inner 
fuel  plasma  transfers  thermal  energy  to  the  outer  propellant  stream  by 
both  convection  and  radiation  heat  transport.  In  this  design,  the  heat 
transfer  occurs  essentially  unimpeded,  but  the  drawback  is  that  the 
propellant  stream  can  physically  mix  with  the  fuel  stream  leading  to  a 
loss  of  nuclear  fuel,  A more  detailed  description  of  the  system  can  be 
found  in  Schneider  and  Thom  (Ref. 6). 

To  avoid  the  mixing  of  the  propellant  and  fuel  streams,  in  the 
nuclear  light  bulb  concept,  the  fissioning  plasma  is  confined  within  a 
transparent  (generally  quartz)  wall.  The  fuel  is  kept  away  from  this 
transparent  wall  by  injecting  a buffer  gas  (argon,  neon,  etc.) 
tangentially  into  the  cavity.  Here,  the  energy  transfer  from  the  plasma 
to  the  propellant  is  by  radiation  heat  transfer  alone.  Major  work  on 
this  concept  was  performed  at  the  United  Aircraft  Research 
Laboratory  (UARL).  Lantham  (Ref,  7)  performed  a series  of  calculations 
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for  the  nuclear  light  bulb  reactor  rocket  engine  under  the  auspices  of 
UARL. 

The  application  of  gas  core  reactors  to  terrestrial  electric  power 
generation  also  had  some  early  investigations.  In  1969,  Gritton  and 
Pinkel  (Ref.  8)  of  the  Rand  corporation  reported  one  such  study  in 
which  they  evaluated  a 4000  MWt  spherical  gaseous  core  system.  The 
central  gas  core  chamber  was  surrounded  by  a moderating/reflecting 
region  and  by  banks  of  energy  conversion  devices.  In  another  study, 
Williams  and  Shelton  of  the  Georgia  Institute  of  Technology  (Ref,9) 
looked  into  the  possibility  of  coupling  gas  core  reactors  with  a 
Magneto  Hydrod3mamic  (MHD)  power  conversion  system  for 
terrestrial  power  generation.  They  studied  the  possibility  of  coupling 
both  the  coaxial  flow  system  and  the  nuclear  light  bulb  system  to  an 
MHD  generator.  They  concluded  that  gas  core  reactors,  in  conjunction 
with  MHD  generators,  are  capable  of  generating  electricity  in  a safe 
and  efficient  manner.  Up  until  then,  the  main  obstacle  preventing 
MHD  generators  from  achieving  their  full  potential  was  the  low 
electrical  conductivity  of  the  working  fluid.  The  study  conducted  at 
the  Georgia  Institute  of  Technology  concluded  that  gas  core  reactors 
are  capable  of  providing  the  high  temperatures  and  the 
non-equilibrium  ionization  necessary  for  high  electrical  conductivities 
and,  hence,  for  high  MHD  efficiencies. 
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Previous  Gas  Core  Reactor  Work  At  The 
University  of  Florida 

As  mentioned  above,  the  major  early  gas  core  reactor  research 
efforts  in  the  United  States  were  closely  related  to  the  NASA's  space 
power  needs  and  the  program  was  essentially  abandoned  in  the  early 
1970s.  But,  at  the  University  of  Florida  (U  of  F),  research  on  gas  core 
reactors  continued  without  much  interruption  during  the  1970s  and 
the  1980s.  During  these  years  many  different  concepts  were  studied 
and  characterized  at  the  U of  F. 

The  work  on  gaseous  core  reactors  started  at  the  U of  F in  the 
early  1970s  when  Schneider  and  Ohanian  (Ref.  10)  proposed  a 
nonsteady  reactor  called  the  Pulsed  Nuclear  Piston  (PNP)  engine.  This 
was  a pulsed  system  operating  on  a thermodynamic  cycle  very  similar 
to  the  internal  combustion  engine  with  the  fissionable  gas  as  the 
working  fluid. 

This  piston  engine  consists  of  a cylindrical  core  enclosed  by  a 
cylinder-piston  arrangement  just  like  in  an  internal  combustion 
engine.  This  is  surrounded  by  a moderating/reflecting  region.  A 
mixture  of  UFe  and  helium  is  used  as  the  fuel  and  the  working  fluid 

(the  helium  is  added  to  enhance  the  thermodynamic  and  heat  transfer 
characteristics  of  the  primary  working  fluid).  The  fuel  is  introduced 
into  the  core,  and  is  then  compressed  by  the  piston  to  form  a 
supercritical  system.  Fission  energy  is  rapidly  released  within  the 
cylinder  thereby  increasing  the  temperature  and  pressure  of  the  fuel 
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considerably.  The  fuel  is  then  exhausted  from  the  cylinder  and  the 
cycle  is  repeated.  Neutronically,  the  core  goes  from  a highly  subcritical 
state  to  a critical  state  and  then  to  a highly  supercritical  state. 

The  energy  released  by  the  fissioning  gas  can  be  extracted  both 
as  mechanical  power  (from  the  motion  of  the  piston)  and  as  heat  from 
the  circulating  gas  (e.g.,  through  a heat  exchanger-turbine 
combination).  Because  of  this  double  mode  of  energy  extraction,  the 
system  efficiency  can  be  as  high  as  50%  (Refill).  An  extensive 
neutronics  and  energetics  analysis  of  this  PNP  concept  was  performed 
by  Dugan  (Ref.  12). 

A pulsed  system  which  is  only  slightly  different  from  the  PNP 
was  also  studied  by  Diaz  et  al.  (Refill)  in  1979.  This  is  the  Pulsed  Gas 
Generator  (PGG)  concept.  The  PGG  is  very  similar  to  the  nuclear 
piston  concept  except  for  the  elimination  of  the  piston  and  making 
the  core  volume  fixed.  The  fuel  (which  is  again  a mixture  of  UFe  and 

helium)  is  cyclically  introduced  into  the  core  just  like  in  a PNP.  The 
pulsing  effect  is  now  achieved  with  the  help  of  rotating  absorber  rods 
placed  in  the  external  moderating  region.  The  PGG  concept  has  the 
advantage  of  mechanical  simplicity  since  it  requires  no  moving  piston. 
However,  the  efficiency  attained  by  the  system  is  lower  than  the  PNP 
due  to  the  fact  that  only  heat  energy  (turbine  power)  is  extracted  from 
the  system  now.  A more  detailed  description  of  the  system  can  be 
found  in  Diaz  et  al.  (Refill). 
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A gaseous  core  system  significantly  different  from  the  pulsed 
systems  (PNP  and  PGG)  was  next  in  line  to  be  analyzed  at  the  U of  F. 
This  is  the  Heterogeneous  Gas  Core  Reactor  (HGCR).  This  novel 
concept  in  gas  core  reactor  design  was  suggested  in  1974  by  Diaz  and 
Dugan  (Ref.  13).  A heterogeneous  core  is  formed  by  arranging  bundles 
of  moderator/ coolant  channels  in  a fissionable  gaseous  core  region. 
Thus  the  HGCR  has  'internal  moderation'  as  opposed  to  the  external 
moderation  in  the  PNP  and  PGG.  The  core  heterogeneity  leads  to 
significant  thermal  hydraulic  and  energetic  advantages  along  with 
improved  technological  feasibility  compared  to  other  gas  core  reactor 
concepts.  Ki  Han  (Ref.  14)  performed  static  neutronic  analysis  and 
some  kinetic  analysis  of  this  system. 

Present  System  Description 

One  of  the  latest  in  the  series  of  gas  core  reactor  concepts  that 
has  been  analyzed  at  the  U of  F is  a bimodal  nuclear  power  system.  A 
schematic  of  the  core  proper  is  shown  in  Fig.  2.1.  Typical  reactor 
dimensions  and  other  key  system  parameters  are  given  in  Table  2. 1 . 
The  development  of  the  system  design  and  characterization  are  being 
done  under  the  auspices  of  the  Innovative  Nuclear  Space  Power 
Institute  (INSPI)  at  the  U of  F. 

As  shown  in  Fig.  2.1,  the  system  consists  of  two  identical  halves. 
Basically,  each  half  of  the  system  consists  of  a central  fissioning  region 
(called  a Burst  Power  Gas  Core  Reactor  or  BPGCR)  surrounded  by  a 
ring  of  mini  cylinders  (called  Pulsed  Gas  Core  Reactors  or  PGCRs) 
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Figure  2-1.  Bimodal  Gas  Core  Reactor  Schematic 


Table  2- 1 . Bimodal  Reactor  Typical  Dimensions  and 
Operating  Conditions. 


PGCR  Data  (Each  Half): 

Number  of  PGCFte 

8 - 12 

Height  (m) 

1.0 

Diameter  (m) 

0.25 

Thermal  Power/PGCR  (MWt) 
BPGCR  Data  (Each  Half): 

3.0 

Diameter  (m) 

1.0 

Height  (m) 

2.0 

Thermal  Power  (MWt) 

200.0 

External  Beryllium  Thickness  (m) 

0.5 

Fuel  (UFe  + He)  Mass  (Kg) 

10.0 

Average  Core  Pressure  (atm) 

20  - 30 

Average  Core  Temperature  (K) 

1500.0 

Helium  Mole  Fraction  in  the  Fuel 

0.93 

U235  Enrichment  in  the  Fuel 

0.85 
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with  a region  of  moderating  beryllium  or  beryllium  oxide  in  between. 
The  PGCRs  are  also  surrounded  by  a reflecting  and  moderating  layer  of 
beryllium  or  beryllium  oxide.  Cooling  channels,  carrying  some  inert 
gas  like  helium,  are  provided  in  the  moderator  region  to  keep  its 
temperature  below  the  melting  point.  Each  half  of  the  core  exhausts 
fuel  gas  into  a central  duct  and  nozzle  configuration  which  then 
directs  the  fuel  gas  into  a disc  MHD  generator.  The  core  fuel  gas 
temperature  is  in  the  range  of  2000-4000  K in  this  concept.  This  high 
temperature  of  the  fuel,  along  with  the  significant  amount  of  fission 
fragment  induced  ionization  taking  place  in  the  MHD  duct,  is 
expected  to  provide  adequate  electrical  conductivity  and  reasonable 
MHD  generator  efficiencies. 

The  PGCRs  are  identical  to  the  pulsed  gas  generators  described 
earlier.  Each  PGCR  consists  of  a cylindrical  core  of  roughly  a meter  in 
height  and  20  to  30  centimeters  in  diameter.  The  core  is  embedded 
in  a beryllium  or  beryllium  oxide  moderating/reflecting  region.  Each 
half  of  the  bimodal  power  system  includes  from  8 to  12  PGCRs,  each  of 
which  is  capable  of  delivering  roughly  3 MWs  of  thermal  power 
(Ref.  15).  The  low  power  mode  of  operation  of  the  bimodal  system  is 
defined,  for  the  purpose  of  this  dissertation,  as  the  case  where  only 
the  PGCRs  are  suppl3dng  power.  The  central  chamber  is  left  empty  or 
filled  with  an  inert  gas.  The  power  from  the  PGCRs  is  expected  to  be 
sufficient  to  meet  the  power  requirements  of  both  the  station  keeping 
mode  and  the  alert  mode  as  mentioned  in  Chapter  1.  Panicker  (Ref.3) 
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has  performed  extensive  static  and  dynamic  neutronic  analysis  of  the 
PGCR  operation  for  this  bimodal  power  system.  Using  Monte  Carlo 
techniques,  he  has  also  studied  the  coupling  effects  among  the  PGCRs 
and  between  the  PGCRs  and  the  BPGCR. 

The  BPGCR  is  a cylindrical  chamber  about  2 meters  in  height 
and  roughly  a meter  in  diameter  (each  half).  This  core  is  surrounded 
by  a beryllium  or  beryllium  oxide  region  which  acts  as  an  external 
moderator  and  reflector.  Each  half  of  the  BPGCR  is  capable  of 
delivering  100s  of  MWt  of  power,  and  during  the  high  power  (or  burst 
power)  mode  of  operation  of  the  bimodal  system,  the  BPGCR  will  be 
supplying  most  of  the  system  power.  The  focus  of  this  dissertation  is 
the  dynamic  analysis  of  the  BPGCR. 

It  should  be  noted  that  since  this  system  is  only  in  a conceptual 
design  state,  variations  of  the  presented  design  are  also  being 
examined.  However,  for  the  purpose  of  this  dissertation,  this  base 
design  is  selected.  The  analysis  that  is  performed  in  this  dissertation 
should  also  be  valid  for  other  designs,  provided  the  concept  is  not 
changed  drastically. 

Summary  of  Previous  Work  Done  in  Related  Areas 
As  mentioned  earlier,  a systematic  study  of  gas  core  reactors  was 
never  undertaken  in  the  United  States,  As  a result  of  this,  plus  the  fact 
that  much  of  the  work  related  to  gas  core  reactors  (space-based 
systems  in  particular)  was  discontinued  in  the  early  1970s,  only  a 
couple  of  studies  on  the  startup  and  stability  of  gas  core  space  nuclear 
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power  systems  are  available.  Although  there  are  many  such  studies  for 
conventional,  terrestrial  nuclear  plants,  they  cannot  be  extended 
directly  to  the  gas  core  space  nuclear  systems  because  of  the  different 
designs  and  quite  different  demands  on  these  systems. 

Gresho  et  al.  (Ref.  16)  and  Birken  (Ref.  17)  deal  with  the  startup 
analysis  of  the  SNAP  reactor  systems  (SNAP  2 and  SNAP  8).  These 
were  small,  solid  core  thermal  reactors  developed  for  meeting  the 
power  requirements  of  various  space  missions.  They  were  capable  of 
generating  only  a few  kilowatts  of  electrical  power.  Since  these 
systems  were  not  meant  for  any  defense  applications,  their  startup 
requirements  were  much  less  severe  than  those  of  the  bimodal  gas 
core  reactor  system  that  is  analyzed  in  this  dissertation.  For  example, 
the  SNAP  2 system  took  about  6 hrs.  in  going  from  a far  subcritical 
state  to  a full  power,  steady  state  operation.  This  long  startup  time 
allowed  the  transition  to  be  made  in  a much  more  controllable  manner 
than  for  the  bimodal  gas  core  system.  For  the  startup  analysis  of  these 
systems,  mostly  analog  computer  techniques  were  used,  and  the  only 
feedback  effects  considered  were  the  fuel  and  grid  plate  temperature 
coefficients.  The  startup  analysis  of  the  gas  core  system  is  expected  to 
be  very  different  from  the  SNAP  system  analysis  because  of  factors  like 
the  gaseous  nature  of  the  fuel,  the  ex-core  circulation  of  the  fuel  gas 
mixture  (which  is  also  the  coolant/ working  fluid),  delayed  neutron 
losses  due  to  the  fuel  circulation,  etc. 
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The  SNAP- 10  reactor,  which  belongs  to  the  same  class  of 
reactors  as  the  SNAP-2,  was  launched  into  space  in  1965.  The  reactor 
was  started  up  while  in  orbit,  and  it  operated  successfully  in  the  space 
environment  for  more  than  a month.  Actual  experimental  data  are 
available  on  the  startup  and  shutdown  of  this  reactor.  A more  detailed 
description  of  the  system  can  be  found  in  Angelo  and  Buden  (Ref. 2). 

In  the  startup  analysis  of  the  nuclear  light  bulb  engine,  the 
authors  (Ref.  18)  did  consider  rapid  startup  situations  (60,  300  and 
600  sec.  power  ramps).  Since  the  objectives  of  the  study  were  to 
estimate  the  thermal  stresses  involved,  the  form  of  fuel  best  suited  for 
the  system,  etc.,  the  neutronic  analysis  of  the  startup  was  performed 
using  only  simplistic  models.  The  analysis  neglected  delayed  neutrons 
completely  and  no  fuel  circulation  effects  were  taken  into 
consideration.  Also,  they  looked  at  only  linear  power  ramps. 

Lantham  et  al.  (Ref.  19)  studied  the  dynamic  response  of  the 
nuclear  light  bulb  engine  to  perturbations  induced  while  the  system  is 
operating  at  full  power.  To  obtain  the  system  response,  they  solved  the 
full  set  of  coupled  thermal,  fluid  dynamics  and  neutron  kinetics 
equations.  Ki  Han  (Ref.  14)  performed  some  stability  analysis  of  the 
high  power  heterogeneous  gas  core  reactor.  His  analysis  was 
preliminary  in  nature  and  he  treated  fuel  circulation  effects  in  a very 
crude  way.  He  also  did  not  include  the  balance-of-plant  in  his  analysis. 


CHAPTER  3 

ANALYSIS  OF  TRANSITION  TO  HIGH  POWER  MODE 

Introduction 

The  purpose  of  the  work  described  in  this  chapter  is  to  analyze 
the  transition  of  the  bimodal  power  system  from  a low  power, 
stationkeeping  mode  of  operation  to  a high  power,  burst  mode.  As 
mentioned  in  Chapter  1,  this  transition  has  to  be  made  very  fast— in 
about  100  seconds— and  it  is  very  important  that  key  system 
parameters  like  core  power,  core  temperature,  core  pressure,  fuel 
mass  flow  rates,  etc.  be  closely  controlled  during  the  transition  to 
make  sure  that  no  maximum  limits  on  any  of  these  variables  are 
exceeded.  The  information  obtained  as  a result  of  this  analysis  can  also 
be  used  in  the  overall  system  design.  For  example,  the  external 
beryllium  (or  beryllium  oxide)  moderator/reflector  has  to  be 
continuously  cooled  to  remove  heat  deposited  in  it  due  to  convective, 
radiative,  and  neutron/gamma  heating  from  the  core  so  as  to  maintain 
this  region  at  the  desired  temperature.  The  heat  deposition  in  the 
moderator  needs  to  be  known  to  establish  the  required  coolant  flow 
rates  through  the  beryllium  to  prevent  it  from  overheating  and 
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melting.  The  temperature  change  of  beryllium  also  introduces 
feedback  reactivity  which  needs  to  be  taken  into  account  when 
determining  the  system's  neutron  multiplication  factor. 

For  the  purpose  of  analyzing  the  transition,  a computer  program 
was  developed.  This  program  can  accept  any  user-specified  reactivity 
variations  as  input,  and  predict  the  behavior  of  the  core  parameters  of 
interest  during  the  transition.  The  program  is  written  in  standard 
FORTRAN  77  and  can  run  either  on  PCs  or  mainframes. 

In  practice,  the  user-specified  input  reactivity  variation  during 
the  transition  can  be  accomplished  in  a number  of  ways.  For  example, 
as  in  a conventional  terrestrial  nuclear  power  plant,  the  fuel  mass 
inside  the  core  can  be  kept  fixed  by  maintaining  constant  and  equal 
fuel  mass  inflow  and  outflow  rates  and  then  with  the  help  of  rotating 
absorber  drums  in  the  external  moderator/reflector,  the  necessary 
reactivity  variation  can  be  brought  about.  Another  way  of  varying  the 
core  reactivity  is  to  vary  the  fuel  mass  inside  the  core.  This  is  possible 
because  of  the  gaseous  nature  of  the  fuel,  plus  the  fact  that  this  is  a 
circulating  fuel  reactor.  The  fuel  mass  variation  inside  the  core  can  be 
achieved  in  a number  of  different  ways.  One  way— which  is  described 
later  in  this  chapter  to  demonstrate  the  working  of  the  computer 
program— is  to  control  the  fuel  mass  inflow  rate  into  the  core  by 
controlling  the  fuel  inlet  pressure  to  the  core.  The  fuel  mass  outflow 
rate  from  the  core  is  then  determined  by  the  conditions  existing 
within  the  core  (temperature,  pressure,  etc.)  and  by  the  downstream 
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conditions.  The  code  can  be  modified,  with  minimum  effort,  to  handle 
any  other  plausible  startup  scenarios. 

System  Description 

The  schematic  diagram  of  the  anal)^ed  reactor  core  is  given  in 
Fig.  2.1.  Typical  dimensions  of  the  core  are  given  in  Table  2.1.  For  the 
purpose  of  this  dissertation,  it  is  assumed  that  the  PGCR's  are  voided 
(unfuelled)  during  the  transition  to  high  power.  A schematic  diagram 
of  the  complete  system,  including  the  balance-of-plant,  is  shown  in 
Fig.  3.1. 

In  Fig.  3.1,  the  two  halves  of  the  BPGCR  core  are  represented  as 
a single  core.  The  output  from  the  core  enters  the  MHD  duct  after 
passing  through  a converging-diverging  nozzle  (not  shown  in  Fig.  3.1). 
The  purpose  of  this  nozzle  is  to  accelerate  the  fuel  exiting  from  the 
core  to  supersonic  velocities  before  entering  the  MHD  generator  (high 
working  fluid  Mach  numbers  are  necessary  for  reasonable  MHD 
generator  efficiencies).  From  the  MHD  duct,  the  fuel  passes  through  a 
supersonic  diffuser  (where  its  velocity  gets  reduced  considerably)  and 
a radiator  heat  exchanger  before  getting  pumped  back  into  the  core.  A 
fuel  reservoir  is  located  before  the  compressor  to  provide  any  excess 
fuel  needed  while  starting  up  the  system,  and  to  act  like  a sink  to 
retain  any  excess  fuel  from  the  loop  during  burst  power  operational 
transients  or  shutdown.  The  heat  exchanger  transfers  heat  to  a 
radiator  for  rejecting  it  to  space. 
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Figure  3-1.  Schematic  of  Complete  BPGCR  Power  System 
Including  Balance-of-Plant. 
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Referring  to  Fig,  2.1,  it  can  be  seen  that  the  reactor  consists  of 
two  identical  halves,  and,  since  the  PGCRs  are  assumed  to  be  shut 
down  during  the  transition  to  high  power,  the  reactor  can  be 
essentially  treated  as  two  large  identical  cylindrical  BPGCR's.  The 
analysis  is  simplified  by  analyzing  the  behavior  of  only  one-half  of  the 
system  during  the  transition.  The  Justification  for  this  simplification  is 
the  assumption  that  the  neutronic  coupling  between  the  two  cores  is 
expected  to  be  small  (this  coupling  actually  depends  largely  on  the 
separation  distance  between  the  two  halves,  i.e.,  on  the  thickness  of 
the  separating  beryllium  or  beryllium  oxide  region).  The  core  power 
and  mass  flow  rates  obtained  from  this  analysis  can  then  be  multiplied 
by  a factor  of  two  to  get  the  corresponding  values  for  the  entire 
system. 

Selection  of  External  Moderator  Thickness 
A 'cavity  reactor'  configuration  is  one  where  the  gaseous  fuel 
region  is  surrounded  by  a moderating/ reflecting  region  and,  unlike  in 
a conventional  solid  core  reactor,  almost  all  of  the  neutron  moderation 
takes  place  in  this  region  external  to  the  fuel  region.  The  BPGCR  is 
such  a cavity  reactor  with  a beryllium  (or  beryllium  oxide)  moderator 
surrounding  the  core.  An  adequate  thickness  of  the  moderating  layer 
is  crucial  for  achieving  criticality.  If  this  thickness  is  too  small,  the 
fission  neutrons  created  within  the  core  do  not  have  a good  chance  of 
being  moderated  in  the  beryllium  region  and  reflected  back  into  the 
core  to  cause  further  fissions.  Most  of  them  will  leak  out  of  the  system. 
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and,  if  the  thickness  is  smcill  enough,  adequate  neutron  thermalization 
also  will  not  occur.  So,  to  maintain  criticality,  either  the  fuel  mass 
inside  the  core  has  to  be  increased,  or  the  core  size  has  to  be 
increased.  At  the  other  extreme,  the  thickness  of  the  beryllium  region 
cannot  be  made  arbitrarily  high  due  to  system  weight  considerations. 
This  is  especially  true  in  space  power  applications  where  one  of  the 
main  goals  of  the  system  design  is  to  keep  the  system  weight  down.  So 
it  is  important  to  find  an  optimum  thickness  for  the  beryllium 
moderator. 

To  find  the  optimum  beryllium  thickness,  the  cylindrical  core  is 
first  mocked  up  in  spherical  geometry  by  conserving  the  total  core 
volume  (see  the  next  section  for  a discussion  on  the  validity  of  this 
procedure).  Then  four  group,  one-dimensional  Sn  transport  theory  keff 

calculations  using  XSDRNPM  (Ref. 20)  are  performed  for  various 
beryllium  thicknesses.  The  obtained  results  are  given  in  Fig.  3.2.  These 
calculations  are  performed  for  an  assumed  core  pressure  of  50 
atmospheres  and  a core  temperature  of  3000  K.  The  fuel  atom 
densities  can  be  estimated  by  assuming  ideal  gas  behavior  for  the  UFe- 

He  mixture.  The  calculated  densities  are  given  in  Table  3.1.  As  can  be 
seen  from  Fig.  3.2,  after  the  beryllium  thickness  exceeds  about  50  cm., 
no  significant  gain  in  system  keff  is  obtained  by  increasing  the 

moderator  thickness  further.  So,  for  the  purpose  of  this  analysis,  a 50 
cm.  beryllium  thickness  is  selected. 


Table  3-1.  Core  Fuel  (UFe-He)  Atom  Densities 


Nuclide 

Atom  Density  (#/cc) 

U235 

7.276  E+18 

U238 

1.283  E+18 

F 

5.135  E+19 

He 

1.137  E+20 

Helium  Mole  Fraction  in  the  Fuel  = 0.93 
Uranium  Enrichment  = 85  atom% 

Average  Core  Pressure  = 50  atm 

Average  Core  Temperature  = 3000  K 


1.20 
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Figure  3-2.  BPGCR  keff  Vs.  Beryllium  Moderator/Reflector  Thickness 
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Spherical  Versus  Cylindrical  Geometry 
Most  of  the  keff  calculations  in  this  dissertation  were  performed 

by  mocking  up  the  actual  cylindrical  geometry  with  an  'equivalent' 
spherical  geometry,  and  then  performing  one-dimensional,  multigroup 
discrete  ordinates  calculation  using  XSDRNPM.  The  'equivalent' 
spherical  geometry  is  formed  by  conserving  the  total  core  volume  and 
the  thickness  of  the  external  beryllium  moderator/reflector.  The 
reason  for  using  a spherical  geometry  instead  of  the  actual  cylindrical 
geometry  is  that  it  is  very  time  consuming  and,  hence,  quite  expensive 
to  perform  a 2-D  or  3-D  cylindrical  geometry  keff  calculation  on  an 

externally-moderated  gas  core  reactor  using  any  standard  neutronics 
code.  One  way  of  getting  around  this  problem  is  by  performing  a 1-D 
(radial  dependence  only)  cylindrical  or  spherical  geometry  calculation. 
When  such  1-D  calculations  are  perfoirned  with  cylindrical 
geometries,  a transverse  buckling  correction  factor  is  generally 
employed  to  account  for  the  perpendicular  (Z-direction)  leakage.  This 
buckling  correction  procedure— which  is  routinely  used  in  solid  core 
reactor  calculations—  is  of  no  value  with  gaseous  core  reactors  because 
of  the  fact  that  the  flux  within  the  core  of  a gaseous  core  reactor  is 
generally  very  flat.  A spherical  geometry  mockup  provides  a better 
one-dimensional  approximation  for  the  cylindrical  gas  core  reactors 
since  the  problem  of  perpendicular  leakage  is  absent  for  this 
configuration. 
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To  check  the  effectiveness  of  the  equivalent  spherical  geometry 
mockup  in  predicating  the  kgff  of  the  actual  cylindrical  core,  a series  of 

calculations  were  performed  using  XSDRNPM  and  MCNP  (Ref.21).  The 
primary  tool  used  for  keff  calculations  in  this  dissertation  is  XSDRNPM. 

It  is  a one  dimensional,  discrete  ordinates  transport  theory  code  from 
Oak  Ridge  National  Laboratory  which  can  perform  multigroup  kgff 

calculations.  MCNP  is  a Monte  Carlo  code  which  can  perform 
continuous  and  discrete  energy  neutron,  photon,  and  coupled 
neutron/photon  transport  through  arbitrary  3-D  geometries.  The 
MCNP  calculations  are  taken  as  the  standard,  and  the  XSDRNPM 
results  are  compared  against  this.  The  XSDRNPM  calculations  were 
performed  using  the  IBM  3090-400  mainframe  computer  at  the 
University  of  Florida  and  the  MCNP  calculations  were  performed  using 
the  CRAY  X-MP/ 48  Supercomputer  at  the  San  Diego  Supercomputing 
Center. 

First,  a set  of  keff  calculations  were  performed  using  XSDRNPM 

and  MCNP  on  a simple,  two  region  spherical  system.  The  system 
consists  of  a spherical  core  of  about  72  cm.  radius  and  a beryllium 
reflector  thickness  of  50  cm.  For  this  set  of  calculations,  both 
XSDRNPM  and  MCNP  treat  the  same  spherical  geometry.  The  results 
of  these  calculation  are  given  in  Table  3.2.  An  identical  geometry 
comparison  was  undertaken  to  eliminate  any  geometry  effects,  and  to 
pinpoint  sources  of  error  arising  from  neutron  cross  section  library 
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Table  3-2.  MCNP  keff  Versus  123-Group,  S4P3  XSDRNPM  kefr 


Geometry 

Calculation  Type 

keff 

% Difference 

Spherical 

Continuous 
Energy  MCNP 

1 

1.1140 

1.01% 

Spherical 

123  Group 
XSDRNPM 

2 

1.1027 

1.  Estimated  Relative  Error  for  MCNP  Result  = 0.7% 

2.  Flux  and  keff  Convergence  Criterion  for 

XSDRNPM  Calculation  = 1.0  E-4 


differences,  including  energy  structure  differences  (see  discussion 
below). 

Calculations  performed  at  the  University  of  Florida,  using 
XSDRNPM,  use  a neutron  cross  section  library  which  comprises  data 
taken  primarily  from  ENDF/B-111  and  ENDF/B-fV  compilations 
whereas  MCNP  uses  mostly  ENDF/B-V  data  maintained  at  the  San 
Diego  Supercomputing  Center. 

The  MCNP  code  is  capable  of  performing  continuous  energy 
neutron  transport  whereas  XSDRNPM  can  perform  only  multigroup 
neutron  transport.  With  the  cross  section  library  maintained  at  the  U 
of  F,  the  finest  group  structure  available  is  123  neutron  groups.  The 
123-group  energy  structure  is  quite  accurate,  and  is  almost  as  good  as 
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a continuous  energy  case.  So,  the  results  given  in  Table  3.2  indicate 
that  the  error  due  to  cross  section  library  difference  is  quite  small 
since  the  1%  or  so  difference  shown  in  this  table  is  within  the 
statistical  fluctuations  of  the  MCNP  results. 

The  next  set  of  calculations  were  performed  to  get  an  estimate 
of  the  error  involved  in  using  a 4-group  energy  structure  compared  to 
a finer  structure,  for  e.g.,  a 26-group  or  a 123-group.  Table  3.3  shows 
the  difference  in  keff  values  obtained  from  a 123-group,  a 26-group, 

and  a 4-group  XSDRNPM  calculations.  The  table  also  gives  the  CPU 
time  taken  for  these  three  cases  on  an  IBM  3090  mainframe.  The 
geometry  and  the  cross  section  library  used  are  the  same  for  all  the 
three  cases.  As  can  be  seen  from  the  table,  the  4-group  calculations 
predicts  a keff  value  which  is  3,1%  higher  than  a 26-group  calculation 

and  the  26-group  calculation  predicts  a keff  which  is  only  0.04% 
higher  than  a 123-group  calculation.  Thus,  it  can  be  seen  that  the 
error  involved  in  using  a 4-group  energy  structure  is  small  enough, 
and  the  computer  time  savings  involved  is  significant  enough  to  justify 
the  use  of  a 4-group  energy  structure. 

Finally,  'equivalent'  spherical  geometry  XSDRNPM  keff 

calculations  (both  123-group  and  4-group)  were  compared  against  a 
cylindrical  geometry  continuous  energy  MCNP  calculation.  The 
'equivalent'  spherical  geometry  was  formed  by  conserving  the  total 
core  volume  and  keeping  the  external  beryllium  thicknesses  the  same. 
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Table  3-3.  XSDRNPM  S4P3  keff  Calculations  - 4-Group  Versus 
26-Group  Versus  123-Group. 


Geometry 

Calculation  TVP^ 

keff 

% Difference 

1 

CPU 

Time  (sec) 

Spherical 

123  Group 

1.1027 

311.48 

Spherical 

26  Group 

1.1031 

0.04% 

21.01 

Spherical 

4 Group 

1.1370 

3.10% 

0.53 

1.  CPU  Time  on  an  IBM  3090-400  Mainframe  for  convergence 
Levels  of  1.0  E-4. 


Table  3-4.  Spherical  Geometry  XSDRNPM  S4P3  keff 
Versus  Cylindrical  Geometry  MCNP  keff. 


Geometry 

Calculation  Type 

keff 

% Difference 

Spherical 

123  Group 
XSDRNPM 

1.1027 

1.2% 

Spherical 

Continuous 
Energy  MCNP 

1.0893 

Cylindrical 

4 Group 
XSDRNPM 

1.1370 

4.4% 

Estimated  Relative  Error  for  MCNP  Result  = 0.7% 
Flux  and  keff  Convergence  Criterion  for 


XSDRNPM  Calculation  =1.0  E-4. 
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The  cylindrical  core  is  2 meters  in  height  and  a meter  in  diameter 
with  a 50  cm.  of  beryllium  surrounding  it.  The  'equivalent'  spherical 
system  has  a radius  of  72. 1 cm.  and  an  external  beryllium  thickness  of 
50  cm.  The  results  of  the  keff  calculations  are  given  in  Table  3.4.  There 

is  a 4.4%  difference  between  the  MCNP  result  and  a 4-group 
XSDRNPM  result.  From  the  previous  two  comparisons  (Tables  3.2  and 
3.3),  we  know  that  for  the  same  geometry,  a 4-group  calculation 
predicts  a keff  value  which  is  about  3%  higher  than  a 123-group 

calculation.  In  addition,  we  also  know  that  the  cross  section  library 
differences  are  responsible  for  small  differences  between  the  keff 

values  predicted  by  XSDRNPM  and  those  predicted  by  MCNP  for  the 
same  geometry.  So,  it  can  be  concluded  that  the  error  due  to  geometry 
difference  (cylinder  vs  sphere)  in  this  particular  case  is  around  2%, 
This  error  should  be  smallest  for  "square"  cylinders  (i.e.,  cylinders 
with  their  height  equal  to  their  diameter),  and  should  increase  as 
cylinders  become  elongated. 

Estimation  of  Gas  Core  Reactor  Reactivity  Coefficients 
Using  the  Equivalent  Spherical  Geometry 

Reactivity  coefficients  are  very  important  neutronics  parameters 
of  a nuclear  reactor,  and,  when  designing  any  new  type  of  reactor, 
these  coefficients— which  are  specific  to  that  system— need  to  be 
determined  carefully.  They  are  used  in  the  stability  analysis  of  the 
system,  and  the  reactor  can  be  stable  or  unstable  at  certain  power 
levels  depending  on  the  magnitude  and  sign  of  these  coefficients. 
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Values  for  these  coefficients  are  also  needed  for  the  startup  analysis  to 
determine  the  total  system  reactivity  during  startup. 

Three  major  reactivity  coefficients  associated  with  the  BPGCR 
are  the  fuel  density  coefficient  (or  equivalently,  the  pressure 
coefficient)  of  reactivity,  the  moderator  temperature,  and  the  fuel 
temperature  coefficients  of  reactivity.  Using  the  selected  equivalent 
spherical  geometry  for  the  gas  core  reactor,  values  for  these  three 
feedback  reactivity  coefficients  were  estimated.  The  results  are  given 
in  Tables  3.5  through  3.8. 

The  fuel  density  coefficient  is  the  most  important  reactivity 
coefficient  for  gas  core  reactors— where  the  fuel  is  in  the  form  of  a 
compressible  fluid.  In  conventional  reactors  (FWRs  or  BWRs),  where 
the  fuel  is  in  solid  form,  this  reactivity  coefficient  is  of  much  less 
importance  because  of  its  much  smaller  magnitude.  The  estimated 
values  of  the  fuel  density  coefficient  of  reactivity  for  the  gas  core 
reactor  are  given  in  Table  3.5.  It  is  a prompt  and  positive  coefficient. 
But  it  should  be  noted  that,  as  a function  of  power  or  temperature,  this 
positive  coefficient  translates  Into  a negative  effect,  i.e.,  if  the  core 
power  or  fuel  temperature  increases,  the  density  of  the  fuel  within  the 
core  is  likely  to  go  down  (since  the  system  is  an  open  flow  system,  an 
increase  in  temperature  need  not  necessarily  imply  a decrease  in  the 
core  fuel  density)  and,  hence,  the  system  reactivity  will  also  go  down. 

The  second  reactivity  coefficient  associated  with  the  fuel  inside 
the  core  is  the  fuel  (Doppler)  temperature  coefficient  of  reactivity. 
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Table  3-5,  Fuel  (UFe-He)  Density  Coefficient  of  Reactivity 


Fuel  Density 
(gm/cc) 

1 

keff 

Fractional  Change  in  keff 
per  Unit  Density  Change 

(Ak/k  per  gm/cc) 

1,0 

0,3752 

+ 450,50 

2,0 

0,5934 

+ 169,50 

4,0 

0,8356 

+ 72,79 

6,0 

0,9668 

+ 40,78 

8,0 

1,0490 

+ 25,50 

10,0 

1,1039 

1.  4-Group  XSDRNPM  S4P3  Calculation  on  a Spherical 

Core  of  Radius  72,1  cm,  with  a Beryllium  Thickness  of  50  cm. 
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Table  3.6  shows  the  system  kefr  as  a function  of  the  fuel  temperature. 
The  fractional  change  in  keff  per  degree  change  in  the  fuel 

temperature  is  also  given  in  this  table.  This  fuel  temperature 
coefficient  of  reactivity  is  a prompt  coefficient  (i.e.,  it  follows  the 
changes  in  the  fuel  temperature  almost  instantaneously)  and.  hence,  is 
a very  important  coefficient  in  conventional  solid  core  reactors  from  a 
safety  standpoint. 


Table  3-6.  Fuel  (UFe-He  Mixture)  Temperature  Coefficient 
of  Reactivity 


Fuel 

Temperature 

(K) 

1 

keff 

Fractional  Change  in  keff 
per  Degree  K 

(Ak/k  per  degree  K) 

500 

1.09867 

-1.766  E-6 

1000 

1.09770 

+4.191  E-7 

1500 

1.09793 

-6.558  E-7 

2000 

1.09757 

1.  123-Group  XSDRNPM  S4P3  Calculation  on  a Spherical  Core 
of  Radius  72.1  cm.  with  a Be  Moderator  Thickness  of  50  cm. 
keff  and  Flux  Convergence  Level=1.0  E-5. 

Increasing  the  fuel  temperature  affects  the  system  reactivity  due 
to  changes  in  the  fuel's  microscopic  cross  section;  this  is  a Doppler 
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effect  in  the  fuel- -where  resonance  interaction  increases  due  to  the 
broadening  of  the  resonances.  In  conventional  reactors,  the  fuel 
enrichment  is  very  low  (3-4%  in  U235)  compared  to  the  80-85% 
enrichment  used  in  the  gas  core  reactors  described  in  this 
dissertation.  With  such  low  enrichment,  the  fuel  temperature 
coefficient  of  reactivity  is  dictated  mainly  by  the  Doppler  effects  in  the 
huge  capture  resonances  of  the  U238  nuclide.  So,  at  all  temperatures, 
the  fuel  temperature  coefficient  is  negative.  In  the  case  of  the  gas  core 
reactor  considered  in  this  dissertation  (where  the  fuel  is  highly 
enriched  in  U235),  the  fuel  temperature  coefficient  is  dictated  by  very 
complex  phenomena,  where  the  Doppler  effect  in  the  capture 
resonances  of  U238  and  U235  compete  with  that  in  the  fission 
resonances  of  U235.  Apart  from  this,  there  is  a neutron  flux  spectral 
effect  also— where  the  increased  resonance  interactions  tend  to 
change  the  neutron  flux  spectral  distribution  within  the  core.  Some 
efforts  were  made  to  pinpoint  the  source  of  the  oscillatory  behavior  of 
the  fuel  temperature  coefficient  of  reactivity  (see  Table  3.6).  But  since 
no  definite  conclusion  could  be  drawn  from  the  limited  study 
undertaken,  this  task  is  left  as  a suggestion  for  future  work. 

A change  in  the  moderator  temperature  affects  the  system 
reactivity  through  two  separate  mechanisms.  The  first  one  is  the 
'spectral  effect’.  Since  almost  all  of  the  neutron  thermalization  is 
taking  place  in  the  beryllium  moderator,  an  increase  in  the  moderator 
temperature  results  in  a harder  thermal  neutron  spectrum  within  the 
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core  (due  to  the  increased  vibrational  energy  of  the  moderator  atoms). 
A harder  thermal  neutron  spectrum  in  turn  implies  a lower  overall 
thermal  neutron  absorption  (fission  plus  radiative  capture)  probability. 
This  part  of  the  beryllium  moderator  temperature  coefficient  of 
reactivity  is  given  in  Table  3.7.  The  values  given  in  this  table  were 
obtained  by  repeating  keff  calculations  in  which  only  the  beryllium 
thermal  cross  section  data  were  changed,  i.e.,  basically  using  beryllium 
cross  sections  with  the  thermal  scattering  kernels  processed  at 
different  temperatures.  Changing  the  thermal  scattering  kernel  of  the 
beryllium  will  alter  the  thermal  neutron  spectrum  within  the  core 
and,  hence,  will  change  the  neutron  multiplication  factor.  All  other 
parameters,  including  the  beryllium  density,  were  kept  constant  for 
these  calculations. 

As  can  be  seen  from  Table  3.7,  the  spectral  component  of  the 
beryllium  temperature  coefficient  of  reactivity  switches  from  a small 
positive  value  at  low  temperatures  to  an  even  smaller,  but  negative, 
value  above  600  K,  and  then  becomes  larger  and  more  negative  with 
increasing  moderator  temperature.  From  a safety  point  of  view,  it  is 
desirable  to  have  this  coefficient  negative  throughout  the  entire 
temperature  range  of  interest.  The  positive  coefficient  at  low 
temperature  is  not  much  of  a concern  because  of  the  fact  that  the 
temperature  range  where  the  coefficient  is  positive  is  quite  small,  and 
at  normal  operating  power  levels  of  the  system,  the  beryllium 
temperature  is  expected  to  be  significantly  higher  than  600  K.  Also,  in 
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practice,  due  to  the  large  thermal  inertia  associated  with  the  beryllium 
(or  beryllium  oxide)  moderator/reflector,  this  reactivity  coefficient  has 
a very  long  time  constant  associated  with  it  and,  hence,  will  contribute 
a negligible  amount  of  feedback  reactivity  during  any  fast  transient. 


Table  3-7.  Spectral  Component  of  Beryllium  Moderator 
Temperature  Coefficient  of  Reactivity. 


Beryllium 
Physical 
Temperature  (K) 

keff^ 

Fractional  Change  in  keff 
per  Degree  K 
(Ak/k  per  degree  K) 

296 

1.08638 

+3.485  E-5 

600 

1.09789 

-2.186  E-6 

800 

1.09741 

-1.549  E-5 

1000 

1.09401 

1.  123  Group  XSDRNPM  S4P3  Calculation  on  a Spherical  Core 
of  Radius  72.1  cm.  with  a Beryllium  Moderator  Thickness  of 
50  cm. 

keff  and  Flux  Convergence  Level  =1.0  E-5. 


The  second  mechanism  through  which  the  moderator 
temperature  affects  the  system  reactivity  is  the  density  effect,  i.e.,  the 
decrease  in  the  physical  density  of  the  moderator  with  increasing 
temperature.  The  beryllium  moderator  density  coefficient  of  reactivity 
is  given  in  Table  3.8.  This  coefficient  has  a negative  effect  with  respect 
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to  temperature  (i.e,,  as  the  physical  temperature  of  the  moderator 
increases,  the  beryllium  density  will  go  down  and,  hence,  the  keff  will 
also  decrease). 


Table  3-8.  Beryllium  Moderator  Density  Coefficient  of 
Reactivity. 


Beryllium 

Physical 

Temperature  (K) 

Beryllium 

Density 

(gm/cc) 

keff^ 

Fractional  Change 
in  keff  per  Degree  K 
(Ak/k  per  degree  K) 

296 

1.8477 

1.08994 

-7.501  E-6 

600 

1.8249 

1.08743 

-9.196  E-6 

800 

1.8063 

1.08543 

-1.230  E-5 

1000 

1.7863 

1.08276 

1.  123-Group  XSDRNPM  S4P3  Calculation  on  a Spherical  Core 
of  Radius  72. 1 cm.  with  a Beryllium  Moderator  Thickness  of 
50  cm. 

Flux  and  keff  Convergence  Level  =1.0  E-5. 


For  obtaining  the  moderator  and  the  fuel  temperature 
coefficients  of  reactivity,  123-group  XSDRNPM  keff  calculations  were 

performed.  For  the  fuel  density  coefficient  of  reactivity  only  4-group 
calculations  were  done.  The  reason  for  selecting  a finer  (and,  hence,  a 
more  accurate)  123-group  calculation  for  the  moderator  and  the  fuel 
temperature  coefficient  determination  is  that  the  magnitudes  of  these 
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coefficients  are  quite  small  (compared  to  the  fuel  density  coefficient  of 
reactivity),  and  if  a crude  energy  structure  were  to  be  used  for 
calculating  them,  it  may  not  be  able  to  pick  up  the  energy  sensitive 
effects  accurately  enough.  Also,  for  the  same  reason,  very  stringent  flux 
and  keff  convergence  criteria  (10-5  or  less)  were  imposed  for  the 

calculation  of  these  two  coefficients.  Compared  to  the  moderator  and 
fuel  temperature  coefficients,  the  density  coefficient  of  reactivity  is 
quite  large  (compare  Table  3.5  with  Tables  3.6,  3.7,  and  3.8)  and, 
hence,  only  4-group  calculations  were  performed  to  obtain  this 
coefficient. 

For  the  transition  to  burst  power  analysis,  only  the  fuel  density 
coefficient  of  reactivity  is  taken  into  account.  The  fuel  temperature 
coefficient  is  neglected  because  of  its  small  magnitude  compared  to 
the  fuel  density  coefficient,  and  the  beryllium  temperature  coefficient 
is  neglected  during  the  transition  because  of  the  long  time  it  takes  for 
the  beryllium  to  undergo  a significant  temperature  change.  The 
change  in  beryllium  temperature  during  the  transition  to  burst  power 
was  calculated  for  various  scenarios  and  is  given  later  in  this  chapter. 

It  was  found  that,  even  for  the  case  where  no  heat  is  removed  from  the 
beryllium  moderator,  the  moderator  temperature  change  during  the 
100  sec.  transition  is  not  significant  enough  to  warrant  the  inclusion  of 
the  beryllium  temperature  coefficient  of  reactivity  in  the  transition.  In 
contrast,  when  the  system  is  operating  at  high  or  burst  power  levels, 
and  full  power  transients  are  analyzed,  the  moderator  temperature 
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coefficient  is  one  of  the  significant— although  slow— feedback  effects. 

So,  in  Chapter  4,  where  the  stability  analysis  of  the  system  is 
described,  this  coefficient  is  included. 

Models  Employed  in  the  Computer  Program 
Modelling  of  the  Core  Geometry 

As  mentioned  earlier,  the  PGCRs  are  assumed  to  be  shutdown 
during  the  transition  to  high  power.  This  makes  it  possible  to  treat  the 
reactor  as  two  identical  halves  (BPGCRs).  The  analysis  is  further 
simplified  by  assuming  a weak  neutronic  coupling  between  the  two 
cores  and,  hence,  only  one  half  of  the  system  is  treated  during  the 
transition.  One  half  of  the  BPGCR  is  shown  in  Fig.  3.3.  The  core  is 
taken  to  be  2 meters  in  height  and  a meter  in  diameter.  It  is 
surrounded  by  a 50  cm.  thick  beryllium  moderator /reflector  (not 
shown  in  the  figure). 

As  shown  in  Fig.  3.3,  the  BPGCR  is  modelled  as  an  'open  flow' 
system  with  no  valves  at  the  core  inlet  or  exit;  in  a real  system,  there 
may  be  some  sort  of  valving  or  orificing  at  the  core  inlet  for  controlling 
the  fuel  flow  into  the  system  and,  if  necessary,  for  shutting  the  system 
down.  Similarly,  there  might  be  a valve  or  orifice  at  the  core  exit  also 
to  control  the  fuel  flow  out  of  the  core.  The  flow  characteristics  will  be 
altered  drastically  if  the  flow  through  the  core  is  controlled  in  this 
manner  and,  hence,  the  fluid  flow  analysis  performed  here  would  have 
to  be  modified  to  take  this  into  account.  A converging  nozzle  is 
provided  at  the  core  inlet.  At  the  core  exit  there  is  a 
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converging-diverging  (C-D)  nozzle.  The  converging  nozzle  at  the  core 
inlet  is  expected  to  provide  an  'inherent'  control  mechanism  while  the 
system  is  operating  at  full  power.  For  example,  if  the  system  power 
level  goes  up  for  some  reason,  the  core  pressure  and  temperature  will 
also  go  up,  and  this  will  cause  a reduction  in  the  fuel  inflow  rate 
(provided  the  flow  at  the  inlet  is  non-choked  and  the  core  inlet 
pressure  is  kept  constant).  The  reverse  action  occurs  if  the  system 
power  drops  for  some  reason.  The  C-D  nozzle  at  the  core  exit  also  will 
provide  a similar  controlling  action.  The  exit  nozzle  is  assumed  to  be 
operating  under  choked  conditions.  The  C-D  nozzle  is  also  needed  at 
the  exit  for  accelerating  the  fuel  to  supersonic  velocities  for  the  MHD 
generator.  The  arrangement  of  the  balance-of-plant  is  shown  in  Fig. 

3.1. 

Thermodynamic  and  Heat  Transfer  Modelling 

The  fuel  (UFe-He  mixture)  is  assumed  to  behave  like  an  ideal  gas 

with  constant  specific  heats  (for  the  purpose  of  this  study,  changes  in 
the  chemical  composition  of  the  fuel  are  neglected).  A lumped 
parameter  approach  is  selected  for  the  thermod)oiamic  and  heat 
transfer  treatment  of  the  system.  For  example,  it  is  assumed  that  an 
average  value  of  the  fuel  temperature  may  be  used  to  describe  the  heat 
transfer  properties  of  the  fuel  inside  the  core  at  any  given  time.  The 
same  is  true  for  thermodynamic  properties  like  pressure,  enthalpy, 
etc..  An  average  value  is  also  assumed  for  the  temperature  and  the  heat 
capacity  of  the  moderator. 
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So,  for  example,  we  can  write  differential  equations  governing 
the  dynamic  behavior  of  the  average  core  temperature  from  a simple 
heat  balance  of  the  form 

Core  Heat  Accumulation  Rate  = Core  Heat  Deposition  Rate 

- Core  Heat  Removal  Rate.  (3-1) 

The  rate  of  heat  accumulation  in  the  fuel  is  given  by 

MHt)  C^  ^ 
dt 

where 

Mf(t)  = instantaneous  fuel  (UFe+He)  mass  inside  the  core  (Kg), 
Cpf  = average  fuel  specific  heat  at  constant  pressure  (J/Kg-K), 
and 

Tf  (t)  = average  instantaneous  core  fuel  temperature  (K). 

The  rate  of  heat  deposition  within  the  core  is  taken  to  be  a constant 
fraction  of  the  instantaneous  core  fission  power,  i.e., 

Qcore  = 0.9  * N (3-2) 

where 

Qcore  = heat  deposition  rate  within  the  core  (watts),  and 
N = instantaneous  core  fission  power  (watts). 

So,  it  is  assumed  that  ninety  percent  of  the  total  fission  energy 
generated  within  the  core  gets  deposited  within  the  core.  The  other 
ten  percent  escapes  from  the  core  in  the  form  of  prompt  and  delayed 
gammas  plus  the  kinetic  energy  of  the  fission  neutrons.  Almost  all  of 
this  energy  is  assumed  to  be  deposited  in  the  moderator/reflector. 
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Three  mechanisms  of  heat  loss  from  the  core  are  considered: 
radiative  transfer  to  the  wall/moderator,  convective  heat  transfer,  and 
the  loss  due  to  the  energy  transported  by  the  flowing  fuel  into  and  out 
of  the  core. 

The  net  radiant  heat  transfer  rate  between  the  core  and  the 
wall /moderator —which  are  taken  to  be  at  uniform  temperature— can  be 
wrritten  as  (see  Ref.  22  for  details) 

9rad  = Ac  O (£  Tf^  - OC  (3-3) 

where 

Qrad  = net  radiant  heat  transfer  rate  between  the  core 
and  the  wall  (watts), 

Ac  = surface  area  of  the  wall/beryllium  moderator  which  is 
exposed  to  the  core  (fuel)  (m2), 
o = Stefan-Boltzmann's  constant  (watts  / m2-K4  ), 

Tf=  average  fuel  temperature  inside  the  core  (K), 

Tw  = average  wall/moderator  temperature  (K), 

8 = average  fuel  emissivity,  and 
a = average  fuel  absorptivity. 

(Implicit  in  the  use  of  the  above  equation  is  the  assumption  that  the 
two  bodies  are  at  uniform  temperature,  i.e.,  radiation  from  each  body 
is  characterized  by  a single  average  temperature.) 
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If  a further  simplifying  assumption  of  a = £ (i.e.,  Kirchhoffs  law) 
is  made,  the  above  equation  becomes 

qrad  = A o e (Tf4  -Tw4).  (3-4) 

(Note  that  in  actual  practice,  radiative  heat  transfer  from  a gas 
presents  a far  more  difficult  problem  than  is  indicated  by  the  above 
equation.  It  is  complicated  by  the  spectral  and  directional  effects 
involved  with  gaseous  absorption  and  emission,  and  also  by  the 
dependence  of  both  a and  £ on  temperature.  But,  for  the  present 
work,  the  simplified  treatment  given  above  is  presumed  to  be 
sufficient.) 

The  film  coefficient  of  heat  transfer,  H,  for  incompressible  pipe 
flow  is  generally  based  on  the  difference  between  the  wall  temperature 
and  mean  stream  temperature.  For  compressible  flow  at  high  speeds, 
it  is  more  appropriate  to  base  H on  the  difference  between  the  actual 
wall  temperature  T^,  and  the  'adiabatic*  wall  temperature.  Taw  (see 

Ref.23  for  details).  Accordingly,  the  total  convective  heat  transfer  rate 
between  the  core  and  the  wall/reflector  can  be  written  as 

Qconv  = H A ( Tw  - Taw  ) (3-5) 

where 

qconv  = net  rate  of  convective  heat  transfer  (watts). 
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H = effective  heat  transfer  coefficient  from  fuel  to  the  wall 
(watts/m2  K),  and 

A = surface  area  of  the  wall/beryllium  moderator  which  is 
exposed  to  the  flowing  fuel  in  the  core  (m2). 

The  adiabatic  wall  temperature  is  defined  by 

1 +R^^M2 
T 2 

where 

T = average  static  temperature  of  the  stream  (K), 

R = recovery  factor  (value  between  0.87  and  0.91  for  subsonic 
pipe  flow), 

M = Mach  number  of  the  flow,  and 
Y = ratio  of  specific  heats  for  the  flowing  gas  (fuel). 

The  film  coefficient  of  heat  transfer,  H,  is  defined  for  turbulent 
flow  (the  flow  through  the  core  will  be  highly  turbulent  even  for 
relatively  small  fuel  velocities  of  the  order  of  20  m/sec)  with  the  help 
of  the  dimensionless  Nusselt  number,  Reynolds  number  and  Prandtl 
number  as  follows: 

Nu  = 0.0364  (Re  * Pr)0  75  (3-6) 

where 

Nu  = Nusselt  number. 

Re  = Reynolds  number,  and 
Pr  = Prandtl  number. 
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(See  Ref. 23  for  details  of  the  above  treatment  of  convective  heat 
transfer.) 

Neglecting  the  changes  in  the  kinetic  and  potential  energies  of 
the  flowing  fluid  across  the  core,  the  energy  loss  due  to  the  energy 
transported  by  the  flowing  fluid  into  and  out  of  the  core  can  be  written 
as 

fltrans  “ r^out  Cp^Tout  " n^ln  Cp^Tin  (3-7) 

where 

niout(t)  = instantaneous  fuel  mass  flow  rate  out  of  the  core 
(Kg/  sec). 

min(t)  = instantaneous  fuel  mass  flow  rate  into  the  core 
(Kg/  sec). 

Tin  = Average  fuel  temperature  at  the  core  inlet  (K),  and 
Tout  = Average  fuel  temperature  at  the  core  exit  (K). 

This  includes  the  energy  loss  associated  with  the  thermal  energy 
changes  across  the  core  plus  the  flow  work  changes. 

Equations  3.2.  3.4,  3.5,  and  3.7  can  be  substituted  in  Eq.  3.1  to 
get  a final  expression  for  the  average  core  fuel  temperature. 

A similar  treatment  is  also  performed  for  finding  the 
instantaneous  average  beryllium  moderator  temperature,  i.e.. 

Be  Heat  accumulation  rate  = Be  Heat  deposition  rate  - 


Be  Heat  removal  rate  (3  - 8 ) 
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The  rate  of  heat  accumulation  in  the  beryllium  moderator  is  given  by 
Mee  CpBe  (dTse  / dt) 

where 

Mbc  = mass  of  beryllium  moderator/ reflector  (Kg), 

CpBe  = average  specific  heat  of  beryllium  (J/Kg-K),  and 
Tbc  = average  beryllium  temperature  (K). 

The  rate  of  heat  deposition  in  beryllium  has  three  components 
to  it.  The  first  one  is  the  heat  deposition  due  to  the  neutrons  and 
gammas  escaping  from  the  core.  This  is  taken  to  be  a constant  fraction 
of  the  instantaneous  core  fission  power.  So, 

QBe  = o.l*N  (3-9) 

where 

QBe  = rate  of  heat  deposition  due  to  neutrons  and  gammas  in 
the  beryllium  (watts),  and 
N = instantaneous  core  fission  power  (watts). 

(That  is,  here  it  is  assumed  that  the  ten  percent  energy  lost  from  the 
core  due  to  escaping  neutrons  and  gammas  is  given  up  fully  in  the 
moderator.  See  also  Equation  3.2.) 

The  second  source  of  heat  deposition  in  the  beryllium 
moderator  is  the  radiative  transfer  from  the  core.  This  part  is  given  by 
Equation  3.4  above.  The  third  source  of  heat  transfer  to  the  beryllium 
is  the  convective  transfer  from  the  flowing  fuel  inside  the  core  to  the 
beryllium  wall.  This  contribution  is  given  by  Equation  3.5. 
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External  Loop  Modelling 

An  ideal  Brayton  cycle  is  taken  to  be  the  thermodynamic  cycle. 
The  T-S  diagram  for  the  cycle  is  shown  in  Fig.  3.4.  The  heat  addition 
process  within  the  core,  and  the  heat  rejection  to  the  heat 
exchanger/radiator  are  assumed  to  take  place  at  constant  pressure. 
The  MHD  generator  is  modelled  as  a conventional  turbine.  The 
performance  of  the  turbine  is  then  dictated  by  selecting  an  enthalpy 
extraction  ratio  and  an  isentropic  efficiency  for  the  turbine  (these  are 
specified  as  constants  to  the  program).  The  values  selected  for  the 
enthalpy  extraction  ratio  and  the  isentropic  efficiency  are  based  on 
some  earlier  analysis  of  the  MHD  generator  performed  by  Gerard 
Welch  at  the  U of  F.  The  enthalpy  extraction  ratio,  together  with  the 
isentropic  efficiency,  can  be  used  to  obtain  values  for  the  stagnation 
pressure  and  temperature  change  across  the  turbine  (MHD  generator). 
The  performance  of  the  compressor  is  handled  by  specifying  an 
isentropic  efficiency— which  is  also  taken  to  be  a constant. 

In  general,  the  analysis  of  the  heat  exchanger/radiator  is  very 
involved  and,  in  the  absence  of  detailed  prior  analysis,  it  cannot  be 
treated  in  a simplistic  way  as  is  done  for  the  MHD  generator  or  the 
compressor  (see  discussion  above).  That  is,  the  behavior  of  the  heat 
exchanger/ radiator  cannot  be  described  by  one  or  two  performance 
indices  as  is  done  for  the  compressor  or  the  MHD  generator. 

To  place  the  treatment  of  the  heat  exchanger/radiator  at  the 
same  level  as  the  analyses  of  the  MHD  generator  and  the  compressor. 
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Figure  3-4.  T-S  Diagram  for  the  BPGCR  Power  System 
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and  also  to  reduce  the  computational  time  for  the  transition  to  burst 
power  analysis,  one  of  two  alternate  simplifying  assumptions  can  be 
made:  either  it  can  be  assumed  that  the  heat  exchanger/radiator 
combination  rejects  heat  at  a constant  rate  or  it  can  be  assumed  that  it 
rejects  a constant  fraction  of  the  heat  generated  within  the  system 
(the  total  heat  includes  the  heat  generated  within  the  core  due  to 
fissions  occurring  there  plus  the  work  done  by  the  compressor  on  the 
system).  For  a gas  thermodynamic  cycle  like  the  Brayton  cycle,  the 
work  done  by  the  compressor  is  quite  significant  (see  discussion 
below).  Note  that  selection  of  either  assumption  implies  that  the 
radiator  has  to  be  of  variable  size  (area).  This  is  because  of  the  fact  that 
during  the  transition  to  burst  power,  the  heat  generated  by  the  system 
(and,  hence,  the  temperature  of  the  working  fluid)  is  continuously 
increasing,  and  if  it  assumed  that  the  heat  rejection  rate  is  constant  or 
is  a constant  fraction  of  the  heat  generated  within  the  system,  then 
the  radiator  size  must  be  changing  continuously.  Such  variable  area 
radiator  concepts  have  been  proposed. 

For  the  purpose  of  this  analysis,  the  second  alternative  is 
selected  since  it  predicts  an  increasing  heat  rejection  rate  during  the 
transition  and,  hence,  is  closer  to  the  actual  situation  than  the  first 
assumption  (i.e.,  in  reality  as  the  system  power  increases,  the  amount 
of  heat  rejected  is  also  expected  to  increase,  and  dramatically  so 
because  of  the  T4  dependence  of  the  radiative  heat  transfer). 
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Bravton  Cycle  vs  Rankine  Cvrle 

For  the  purpose  of  this  dissertation,  the  thermod5mamic  cycle  is 
taken  to  be  an  ideal  Brayton  cycle.  The  primary  motivation  for 
emplo)dng  a gaseous  cycle  like  the  Brayton  cycle  was  that  the 
fuel/working  fluid  considered  initially  was  a UFe-He  mixture,  and  at 

typical  operating  temperatures  and  pressures  existing  within  the 
system,  this  mixture  always  remained  in  a gaseous  form.  Also,  initially, 
the  possibility  of  a pulsed  system  with  acoustic  amplification  was 
considered  for  the  burst  power  mode  of  operation.  It  was  conjectured 
that  under  these  conditions,  a Bra)d:on  cycle  might  prove  to  be 
satisfactory.  INSPI  has  also  looked  into  the  possibility  of  using  a 
fuel/working  fluid  mixture  consisting  of  UF4  and  helium.  The  UF4  can 

be  condensed  out  of  the  mixture  at  the  temperatures  and  pressures 
existing  within  the  system  and,  hence,  the  possibility  of  a hybrid 
Rankine /Brayton  thermod)mamic  cycle  was  also  considered.  Currently, 
INSPl  is  investigating  fuel/working  fluid  mixtures  of  UF4  and  metal 

fluorides  in  which  the  entire  fuel/working  fluid  mixture  is  condensible 
under  expected  operating  conditions.  Such  a fuel/working  fluid 
mixture  means  that  the  system  can  operate  on  a true  Rankine  cycle. 

The  main  drawback  of  the  Brayton  cycle  is  the  compressor  work 
needed  to  compress  the  gas.  This  is  significantly  higher  for  the 
Brayton  cycle  than  for  a Rankine  cycle.  In  the  latter,  since  only  an 
almost  incompressible  liquid  is  pumped,  the  typical  hackwork  ratio 
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(which  is  defined  as  the  ratio  of  the  pump  or  compressor  work  to  the 
turbine  work)  is  only  a few  percent  or  less  compared  to  50-80%  for  a 
Brayton  cycle  (for  an  ideal  Bra)d;on  cycle  it  is  of  the  order  of  50%  but 
in  actual  systems,  it  can  be  as  high  as  80%). 

A second  drawback  of  the  Brayton  cycle,  which  is  very  significant 
for  space  power  systems,  is  the  temperature  at  which  heat  rejection 
occurs.  In  a Rankine  cycle,  the  heat  rejection  can  take  place  at 
constant  temperature  (condensing  vapor)  whereas  in  a Brayton  cycle  it 
occurs  at  a varying  temperature.  If  we  were  to  determine  an  "effective" 
constant  temperature  for  heat  rejection  for  a Brayton  cycle,  it  turns 
out  to  be  much  lower  than  the  constant  temperature  at  which  the 
Rankine  cycle  rejects  heat  (provided  we  place  comparable  thermal 
limits  on  the  Brayton  and  Rankine  cycles).  This  factor  can  dramatically 
increase  the  size  of  the  radiator. 

If  the  frictional  losses  in  the  system  are  neglected,  and  if  the 
efficiencies  of  a Brayton  cycle  and  a Rankine  cycle  are  compared  (for 
roughly  the  same  operating  conditions),  the  Rankine  cycle  can  actually 
turn  out  to  be  poorer  than  a Brayton  cycle  (22-26%  for  a Brayton  cycle 
compared  to  roughly  22%  for  a Rankine  cycle)  (numbers  are  based  on 
some  earlier  work  done  by  Gerard  Welch  at  the  U of  F)).  However,  we 
cannot  neglect  frictional  loses,  and  making  up  for  these  in  a Rankine 
cycle  costs  almost  nothing;  the  cost  in  a Brayton  cycle  can  be  very  high 
because  of  the  high  hackwork  ratio.  With  the  frictional  loses  included, 
the  Brayton  cycle  efficiency  turns  out  to  be  lower  than  that  for  the 
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Rankine  cycle.  This  lower  efficiency  also  translates  into  a bigger 
radiator  mass  (since  a lower  efficiency  implies  more  heat  rejection  for 
the  same  electrical  power  output). 

Fluid  Flow  Modelling 

As  mentioned  earlier,  the  fuel  (UFe-He  mixture)  is  assumed  to 

behave  like  an  ideal  gas  with  constant  specific  heats.  The  fuel  is  also 
assumed  to  maintain  its  integrity  throughout  the  transition  to  burst 
power  (i.e.,  no  dissociation  of  the  UF©  is  taken  into  account). 

Even  after  the  commencement  of  the  transition,  it  takes  some 
time  before  the  core  temperature  and  pressure  builds  up  to  levels 
where  significant  dissociation  of  UFe  can  occur.  Significant  UFe 

dissociation  occurs  only  at  temperatures  above  1600  K at  atmospheric 
pressure.  With  increasing  gas  pressure,  dissociation  becomes 
significant  only  at  higher  temperatures.  For  e.g.,  with  a UFe  partial 

pressure  of  about  10  atm.,  the  fuel  can  exist  without  much  UFe 

dissociation  up  to  about  2200  K.  The  fuel  temperature  profile  during 
the  transition  is  given  in  Fig.  3.21  later  in  this  Chapter.  From  this 
figure,  it  can  be  seen  that,  at  least  for  the  first  50  sec.  of  the  total  100 
sec.  transition  time,  the  fuel  temperature  remains  below  the  2200  K 
value.  (UFe  may  be  able  to  exist  without  significant  dissociation  at 

much  higher  temperatures--up  to  about  3000  K,  provided  fluorinating 
agents  like  CIF3  or  BiFa  are  added  to  the  fuel  (Ref  24)). 
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If  the  fuel  is  in  the  form  of  UF4.  the  dissociation  is  not  at  all  a 

concern  even  at  temperatures  as  high  as  3000-4000  K.  Neutronically, 
it  doesn't  matter  whether  the  fuel  is  in  the  form  of  UF4  or  UFe.  The 

extra  fluorines  have  no  significant  neutronic  effect.  So,  the  neutronic 
analysis  done  here  remains  valid  whether  UFe  or  UF4  is  used  as  the 

fuel. 

Core  hydrodynamics  is  represented  by  fully  developed, 
one-dimensional,  turbulent  flow  (as  mentioned  earlier,  the  fuel  flow 
through  the  core  will  be  highly  turbulent  even  for  relatively  low  fuel 
velocities  of  the  order  of  20-50  m/sec).  Frictional  effects  are 
neglected.  Isentropic  flow  is  assumed  for  the  flow  through  the  inlet 
and  exit  nozzles  so  that  stagnation  temperatures  and  pressures  of  the 
fuel  will  adequately  describe  the  flow  of  the  fluid  through  the  nozzles. 
For  example,  the  maximum  fuel  flow  rate  through  the  exit  nozzle  is 
given  by  the  relation 

u,  _ 0.6847  Po  A 

i^max TTo 

R (To)^/2  (3-10) 

where 

nimax  = maximum  fuel  mass  flow  rate  (Kg/sec), 

A = exit  nozzle  throat  area  (m2) , 

Po  = core  stagnation  pressure  (Pa),  and 

To  = core  stagnation  temperature  (K). 
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(This  maximum  flow  occurs  under  choked  flow  conditions,  i.e.,  when 
the  Mach  number  at  the  throat  is  unity).  The  exhaust  nozzle  is 
assumed  to  be  under  choked  flow  conditions  always  (flow  through  the 
converging-diverging  nozzle  at  the  core  exit  has  to  be  choked  to  obtain 
supersonic  flow  going  into  the  MHD  generator).  The  converging  nozzle 
at  the  core  inlet  is  operated  under  non-choked  conditions  for  reasons 
given  earlier  in  the  geometry  modelling  section. 

Neutron  Kinetics  Modelling 

It  was  mentioned  earlier  that  it  is  assumed  that  negligible 
decomposition  of  UFe  takes  place  during  the  transition  to  burst  power. 

Apart  from  this,  it  is  also  assumed  that  the  fuel  composition  remains 
constant  while  going  through  the  external  loop.  This  implies  that  the 
UFe  mole  fraction  of  about  93%  and  the  85%  uranium  enrichment 

remains  the  same.  It  is  also  assumed  that  there  is  no  significant 
buildup  of  strongly  neutron  absorbing  fission  products.  These 
assumptions  should  not  result  in  too  much  error  since  the  fuel  bumup 
and  fission  product  buildup  are  going  to  be  very  small  during  the  short 
transition  time. 

To  model  exactly  the  time  behavior  of  the  neutron  population 
within  a reactor,  use  of  the  time-dependent  Boltzmann's  neutron 
transport  equation  has  to  be  made.  This  is  all  the  more  true  for  the 
type  of  gas  core  reactor  that  is  anal}^ed  in  this  dissertation  because  of 
the  fact  that  within  the  core  proper,  almost  no  neutron  moderation 
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(scattering)  is  taking  place,  and  the  core  is  only  of  the  order  of  a mean 
free  path  or  less  in  size.  The  only  significant  interaction  taking  place 
within  the  core  is  the  absorption  of  thermal  neutrons  returning  to  the 
core  from  the  reflector/ moderator.  But,  from  a practical  standpoint, 
implementing  this  model  can  be  very  tedious  and  computationally 
intensive,  and  is  attempted  only  in  major  time  dependent  transport 
theory  computer  codes. 

Some  previous  studies  conducted  at  the  U of  F indicate  that,  for 
gas  core  reactor  systems  which  posses  a significant  degree  of 
grayness,  multigroup  diffusion  theory  is  generally  adequate  for  static 
calculations  provided  good  fast  and  thermal  group  constants  are  used 
(Ref.  25)  (grayness  for  an  externally  moderated  gas  core  reactor  is 
defined  as  the  ratio  of  the  thermal  neutron  current  into  the  core  from 
the  moderating/reflecting  region  to  the  value  of  the  thermal  neutron 
flux  at  the  core-moderator/reflector  boundary).  For  small  transients 
(i.e.,  small  reactivity  variations),  we  can  assume  that  time-dependent, 
multigroup  diffusion  theory  is  adequate  to  describe  the  dynamic 
behavior  of  the  gas  core  reactor.  Sometimes,  even  this  model  is  too 
detailed  for  practical  implementation  in  reactor  kinetics  analysis  due 
to  excessive  computational  requirements. 

For  this  reason,  in  this  dissertation  (as  is  the  common  practice 
for  many  preliminary  level  reactor  kinetics  studies),  it  is  assumed  that 
the  neutron  dynamic  characteristics  of  the  gas  core  reactor  can  be 
adequately  described  by  the  average  neutron  density,  'n',  i.e.,  the 
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spatial  dependence  of  the  neutron  population  is  neglected.  Often,  in 
neutron  kinetics  studies,  'n'  is  assumed  to  represent  only  the  neutrons 
in  the  energy  range  where  most  of  the  fissions  occur— for  example, 
thermal  neutrons  in  a thermal  reactor.  This  type  of  analysis  is  by  no 
means  exact  but,  for  preliminary  scoping  purposes,  this  is  adequate.  It 
is  capable,  for  example,  of  providing  the  time  dependence  of  the  gross 
reactor  power  level  during  transients. 

With  the  above  assumption,  the  d5mamic  behavior  of  the  neutron 
density  within  the  core  can  be  represented  by 


dn(t) 

dt 


p(t)  - (3eff 

r 


n(t)  + X 

i 


where 

n(t)  = core  neutron  density  at  time  t (#/cm3), 
p(t)  = time  dependent  system  reactivity. 


Peff  = effective  delayed  neutron  fraction. 


(3.11) 


1*  = mean  neutron  generation  time  (sec), 

Ci(t)  = effective  delayed  neutron  precursor  density  for  the  ith 

delayed  neutron  group  (#/cm3),  and 
A.1  = decay  constant  for  the  ith  delayed  neutron  precursor 

group  (sec-i). 

To  complete  the  description  of  the  neutron  dynamics  of  the 
reactor,  an  equation  for  the  ith  delayed  neutron  precursor  density,  Ci, 

also  has  to  be  provided.  For  the  gas  core  reactor  analyzed  here,  this 
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precursor  density  equation  is  quite  different  from  that  for  a 
conventional,  solid  core  reactor.  Since  this  is  a circulating  fuel  reactor, 
not  all  of  the  delayed  neutron  precursors  will  decay  within  the  core 
and  thereby  contribute  to  the  total  core  neutron  population.  Some  of 
the  precursors  are  going  to  be  swept  out  of  the  core  (along  with  the 
fuel)  before  they  can  decay.  Now  two  things  can  happen  to  these 
precursors  which  have  left  the  core  without  decaying;  they  can  either 
decay  outside  the  core  (so  these  neutrons  are  lost  forever,  and  will  not 
contribute  to  the  neutron  population  within  the  core)  or  they  can 
reenter  the  core— after  going  through  the  external  loop  without 
decaying.  In  either  case,  the  effective  fraction  of  delayed  neutron  for 
this  type  of  reactor  is  going  to  be  less  than  that  for  a stationary  fuel 
reactor.  The  precursor  density  equation  taking  into  account  the  fuel 
circulation  is 

dCj(t)  ^ Peff  n(t)  _ ^ , _ ci(t)  ci(t-xi)  exp  (-Xi  ^j) 

dt  1*  ^ ^ (3-12) 

for  i = 1,2 6. 

where 

Xc  = mean  core  residence  time  for  the  fuel  (sec),  and 

Ti  = mean  external  loop  circulation  time  for  the  fuel  (sec) 

(All  other  variables  are  as  defined  in  Equation  3.11.) 

The  third  term  on  the  right  hand  side  of  the  above  equation  accounts 
for  the  loss  of  the  delayed  neutron  precursors  from  the  core  due  to 
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flow,  and  the  last  term  on  the  right  hand  side  takes  care  of  the  reentry 
of  the  precursors  at  the  core  inlet  after  going  through  the  external 
loop. 

So,  the  effect  of  the  fuel  flow  outside  the  core  is  to  make  the 
reactivity  needed  (and,  hence,  the  fuel  mass  needed  inside  the  core) 
for  steady  state  operation  larger  than  the  value  needed  if  it  were  a 
non-flow  system. 

In  Equation  3.11,  'n'  represents  the  instantaneous  neutron 
density  within  the  core.  But  it  can  represent  any  integral  or  volume 
averaged  property  of  the  reactor  that  is  proportional  to  the 
instantaneous  neutron  density.  For  example,  it  can  be  redefined  so 
that  it  represents  the  total  number  of  neutrons  within  the  core  at  any 
given  instant  or  the  fission  rate  or  the  total  core  power  at  any  given 
instant,  etc..  For  convenience,  the  dependent  variable  in  Equation  3.11 
is  defined  to  be  the  instantaneous  power,  N,  of  the  reactor.  Now  the 
precursor  density,  Ci,  used  in  Equations  3.11  and  3.12  has  to  be 

replaced  by  a new  quantity  Ci  defined  as  follows: 

Ci  = WfVXfCiv  (3-13) 

where 

Cl  = effective  delayed  neutron  precursor  density  of  the  ith 
delayed  neutron  group  (#/cm3), 

Wf  = usable  energy  produced  per  fission  (MeV/fission), 
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Zf  = average  thermal  macroscopic  fission  cross  section  (cm- 1 ) , 

V = total  core  volume  (cm3),  and 

V = neutron  velocity  (cm/sec). 

As  a final  simplification,  instead  of  the  usual  six  delayed  neutron 
groups,  only  one  effective  group  of  delayed  neutrons  is  considered,  i.e., 
the  summation  term  in  Equation  3.11  is  replaced  by  one  term  which 
represent  an  "effective"  delayed  neutron  group.  Also,  Equation  3.12— 
which  is  a set  of  six  equations,  is  now  replaced  by  one  equation  for  the 
"effective"  delayed  neutron  group.  This  simplification  is  also  done  so  as 
to  reduce  the  computational  time.  The  computer  program  developed 
can  equally  well  handle  all  six  delayed  groups  (with  a corresponding 
increase  in  the  computational  time). 

There  are  different  ways  of  defining  an  "effective"  one  delayed 
group  of  neutrons.  The  recipe  used  in  this  dissertation  for  obtaining 
the  value  for  J3,  and  the  corresponding  value  for  X is  as  follows; 

.6  = I^i  (3-14) 

and 


(I  Pi/  Xi) 

i 


(3-15) 


for  i = 1,2 6 
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where 

= fraction  of  the  ith  delayed  neutron  group  emitted  per 
fission,  and 

= decay  constant  for  the  delayed  neutron  precursor  for  the 
ith  delayed  neutron  group. 

The  neutron  generation  time,  T,  used  in  Equations  3.11  and  3.12 
is  an  important  factor,  and  is  defined  as  the  average  time  between  the 
birth  of  a neutron  and  its  subsequent  absorption  inducing  fission.  It  is 
related  to  the  prompt  neutron  lifetime,  T of  the  reactor  by  the 
relation 

1*  = 1 / keff 

where  'keff  is  the  effective  neutron  multiplication  factor  of  the  reactor. 

The  prompt  neutron  lifetime  (which  is  defined  as  the  average  time 
between  the  birth  of  a neutron  and  its  subsequent  removal  from  the 
system— either  by  absorption  or  by  leakage)  is  a very  difficult 
parameter  to  estimate  anal5dically  since  it  depends  in  a complex 
manner  on  a number  of  system  variables  like  geometry,  fuel 
enrichment,  reflector  thickness,  etc.  So,  for  the  purpose  of  this  work, 
T for  the  system  is  determined  by  independent  Monte  Carlo 
calculations,  and  the  value  thus  obtained  is  used  as  an  input  to  the 
program. 
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Programming  Aspects 

The  computer  program,  DYNAM,  developed  as  a part  of  this 
work  for  analyzing  the  transients,  is  written  in  standard  FORTRAN  77. 
It  can  run  either  on  a PC  or  a mainframe  (the  program  calls  for  certain 
timing  subroutines  which  are  only  available  on  IBM  mainframes.  These 
call  statements  need  to  be  deleted  from  the  program  before  tr3dng  to 
run  it  on  a PC  or  on  a non-IBM  mainframe). 

Simple  finite  difference  techniques  are  used  for  solving  most  of 
the  relevant  equations  like  those  for  core  temperature,  beryllium 
temperature,  core  fuel  mass,  etc..  The  coupled  neutron  kinetics 
equations  are  solved  using  a fourth  order  Runge-Kutta  method  (this 
method  is  selected  because  it  is  quite  easy  to  implement,  plus  it  has 
good  numerical  stability  characteristics). 

The  only  difficult  term  to  handle  in  the  neutron  kinetics 
equation  is  the  precursor  reintroduction  term,  i.e.,  the  last  term  on 
the  right  hand  side  of  Equation  3.12— which  represents  the  precursors 
which  are  reintroduced  into  the  core  after  going  through  the  external 
loop.  This  term  is  complicated  because,  at  any  given  instant,  it 
depends  on  the  precursor  density  which  had  existed  within  the  core 
at  a time  Ti  (the  mean  external  loop  residence  time)  before. 

When  the  core  power  is  varying  rather  slowly,  a reasonable 
approximation  would  to  be  take 


69 


Cl  (t-Ti)  = Cl  (t) 

This  is  certainly  not  true  during  the  transition  to  burst  power— where 
the  power  goes  from  almost  zero  to  a few  hundred  megawatts  in  about 
100  sec.  So,  in  DYNAM,  the  precursor  reintroduction  term  is  handled 
numerically  as  described  below. 

The  precursor  density  at  each  time  step  is  stored  in  a stack, 
and,  at  any  instant,  the  value  of  'C  which  existed  at  a time  li  before  is 

used  to  estimate  the  precursor  reintroduction  term.  Once  a value 
C(t-Xi)  is  used  during  a time  step,  it  can  be  discarded  from  the  stack  of 

stored  values  of  'C,  That  way,  at  any  given  instant,  the  number  of  values 
of  'C  that  need  to  be  stored  is  only 
Xi  / step 

where  'step'  is  the  (constant)  time  step  size  used  in  the  numerical 
computation. 

Sometimes  this  can  create  memory  limitation  problems  while 
running  on  a PC.  In  this  case,  the  precursor  density  can  be  averaged 
over  a fixed  number  of  time  steps  before  storing  it.  The  number  of 
times  steps  over  which  the  precursor  density  is  averaged  depends  on 
the  transient.  If  it  is  a very  fast  transient,  one  might  want  to  average 
the  density  over  a small  number  of  time  steps,  whereas,  for  slower 
transients,  the  density  can  be  averaged  over  a large  number  of  time 
steps  without  incurring  serious  error. 
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Testing  of  the  Computer  Program 
Before  the  developed  computer  program,  DYNAM,  can  be  used 
for  analyzing  the  transition  to  burst  power,  it  has  to  be  tested  or 
benchmarked.  The  complete  program- -which  includes  the  neutron 
kinetics,  the  thermodynamic,  the  heat  transfer,  and  the  fluid  flow 
aspects  of  the  gas  core  reactor,  cannot  be  benchmarked  against  any 
standard  codes  since  there  are  no  existing  codes  which  are  capable  of 
analyzing  transients  for  such  a novel  type  of  reactor  (novelty  includes 
the  combination  of  gaseous  fuel  in  a closed  circulating  loop,  highly 
compressible  fuel  which  also  serves  as  the  working  fluid,  etc.).  So,  it 
was  decided  to  test  individual  sections  of  DYNAM  against  some 
standard  codes  or  against  analytical  results— where  possible. 

One  of  the  major  subroutines  in  DYNAM  is  the  one  which 
estimates  the  effective  neutron  multiplication  factor  of  the  system  (or 
the  system  reactivity).  This  subroutine  accepts  the  reactor  dimensions 
(core  radius  and  beryllium  thickness)  and  the  fuel  atom  densities  as 
input,  and  it  estimates  the  system  keff.  The  system  keff  needs  to  be 

known  at  every  instant  during  the  transition  to  burst  power. 

Since  the  system  keir  has  to  be  estimated  a large  number  of 

times  during  the  transition  analysis,  this  subroutine  was  written  giving 
speed  of  calculation  priority  over  accuracy.  The  subroutine  mocks  up 
the  cylindrical  core  with  a spherical  geometry,  and  then  performs  a 
two-region  (core  and  moderator/reflector),  two-group  (fast  and 
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thermal)  diffusion  theory  kgff  calculation.  The  microscopic  group 

constants  needed  for  the  keff  calculation  are  provided  as  input  to  the 
subroutine.  These  group  constants  are  obtained  from  independent, 
two-group  transport  theory  calculations  using  XSDRNPM.  These 
microscopic  constants  are  found  to  be  not  very  sensitive  to  small 
changes  in  fuel  atom  densities  within  the  core  or  to  changes  in  the 
core  dimensions.  But,  if  either  of  the  two  (geometry  or  fuel  atom 
densities)  is  changed  significantly,  the  group  constants  need  to  be 
reevaluated  using  XSDRNPM.  The  effects  of  temperature  on  the 
microscopic  constants  are  handled  via  feedback  coefficients.  These 
coefficients  are  given  in  Tables  3.6  and  3.7. 

The  system  keir  values  calculated  by  the  DYNAM  subroutine  are 

compared  against  those  estimated  by  XSDRNPM  (for  the  same 
geometry  and  fuel  atom  densities).  Comparisons  are  made  against 
two-group  diffusion  theory  and  two-group  transport  theory  results 
from  XSDFdMPM.  The  results  of  the  comparisons  are  shown  in  Table 
3.9.  The  comparisons  are  done  for  various  fuel  atom  densities  within 
the  core. 

As  can  be  seen  from  the  table,  in  all  cases,  the  error  involved  is 
less  than  a percent  and,  hence,  is  quite  acceptable— especially 
considering  the  fact  that  the  subroutine  is  quite  fast.  It  takes  less  than 
0,01  sec.  to  estimate  the  keff  of  the  system  to  a convergence  level  of 

better  than  10-4  on  an  IBM  3090  mainframe.  Compared  to  this,  a 
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Table  3-9.  Comparison  of  keff  Values  Obtained  Using  the  kecr 
Subroutine  of  DYNAM  and  the  XSDRNPM  Code  for 
Different  Fuel  Density  Factors. 


1 

Core 

Fuel 

Density 

Factor 

keff 

Percentage  Difference 
Between  DYNAM  keff 

and 

2 

DYNAM 

3 

2-Group 

XSDRNPM 

(Diffusion 

Theory 

Option) 

4 

2 -Group 
XSDRNPM 
(S4P^ 

Two  Group 
XSDRNPM 
(Diff.  Theory) 

Two  Group 
XSDRNPM 
[S^P^ 

0.2 

0.6990 

0.7025 

0.7000 

-0.498 

-0.143 

0.5 

1.0356 

1.0411 

1.0384 

-0.528 

-0.269 

0.8 

1.1760 

1.1789 

1.1758 

-0.246 

+0.017 

1.0 

1.2310 

1.2337 

1.2302 

-0.219 

+0.065 

1 . A fuel  density  factor  of  unity  corresponds  to  the  following  atom 
densities  (in  atoms/barn-cm); 


U235  = 1.458  E-5 

U238  = 2.548  E-6 
F = 1.027  E-4 

He  = 2.275  E-4 

2.  Two  group,  two  region,  diffusion  theory  calculation  on  a 
spherical  core  with  72. 1 cm.  core  radius  and  50  cm.  beryllium 
moderator/reflector  thickness. 

3.  Two  group,  two  region  diffusion  theory  calculation  on  a 
spherical  geometry  of  the  same  dimensions  as  given  above. 

4.  Two  group,  two  region  transport  theory  calculation  (S4P3) 
on  a spherical  system  of  the  same  dimensions  as  given  above. 
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two-group,  two-region  XSDRNPM  keff  calculation  on  a similar  system 
takes  about  0.5  sec. 

The  fact  that  the  DYNAM  results  are  closer  to  the  two  group 
transport  theory  results  is  purely  coincidental.  In  general,  it  was  found 
that  diffusion  theory  predicts  a keir  which  is  higher  than  the  keff 

prediction  of  transport  theory.  So,  the  XSDRNPM  two  group  transport 
theory  results  are  lower  than  the  XSDRNPM  diffusion  theory  results.  It 
so  happens  that  DYNAM  predicted  keff  values  are  also  slightly  lower 

than  the  XSDRNPM  two  group  diffusion  theory  results  and,  thus, 
actually  a little  closer  to  the  two  group  transport  theory  results. 

To  complete  the  benchmarking  of  the  keff  subroutine  of  DYNAM, 

the  fast  and  thermal  group  fluxes  obtained  using  the  keff  subroutine  are 

also  compared  against  those  obtained  using  two  group  transport  theory 
calculations  (XSDRNPM).  The  results  of  this  comparison  are  given  in 
Figs.  3,5  and  3,6.  As  can  be  seen  from  the  figures,  except  for  the  case 
of  fast  flux  within  the  core,  DYNAM  predictions  of  the  thermal  and  fast 
fluxes  are  in  agreement  with  the  predictions  of  XSDRNPM.  DYNAM 
underpredicts  the  fast  flux  within  the  core  by  roughly  10%.  The  fluxes 
shown  in  these  figures  (group  1 and  group  2)  are  normalized  such  that 
the  total  fission  neutron  source  in  the  core  is  one  neutron  per  second 
(it  should  be  mentioned  here  that  the  agreement  between  flux  shapes 
predicted  by  the  diffusion  theory  and  the  transport  theory  is  not 
always  this  good.  For  certain  fuel  density  regimes,  the  agreement 
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Figure  3-5.  Comparison  of  Radial  Flux  (Fast)  Profile  Within  the  BPGCR 
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Figure  3-6.  Comparison  of  Radial  Flux  (Thermal)  Profile  Within  the 
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between  Sn  and  diffusion  theory  fluxes  can  be  very  poor  even  though 
the  keff  values  are  in  reasonable  agreement). 

A second  major  subroutine  in  DYNAM  is  the  one  which  solves 
the  coupled  point  reactor  kinetics  equations.  This  subroutine  uses  a 
fourth  order  Runge-Kutta  method  to  solve  the  coupled  neutron  and 
precursor  density  equations.  Two  separate  methods  were  used  to 
check  the  validity  of  this  subroutine.  First,  the  effectiveness  of  this 
subroutine  in  predicting  the  neutron  density  behavior  after  step 
insertions  of  reactivity  was  compared  against  analytical  results.  For  two 
different  step  reactivity  insertions,  the  results  obtained  using  the 
developed  subroutine  are  compared  against  the  "exact"  analytical 
results  (see  Ref.  26),  and  also  against  the  'prompt  jump’  approximation 
results  (the  prompt  jump  approximation  assumes  that  immediately 
following  a step  change  in  reactivity,  the  delayed  neutron  precursor 
concentration  remains  unchanged  over  the  time  of  the  sudden  rise  or 
drop  in  flux).  The  "exact"  analytical  result  is  for  one  delayed  neutron 
group:  it  gives  good  results  when  the  reactivity  step  is  no  greater  than 
about  0.3  J3  (see  Ref.26  for  the  limitations  of  this  "exact"  anal}d;ical 
expression).  Figures  3.7  and  3.8  gives  the  results  of  this  comparison. 
Again,  it  can  be  seen  that  the  results  from  the  developed  subroutine 
agree  quite  well  with  analytical  predictions. 

Next  the  developed  subroutine  is  compared  against  the  standard 
point  reactor  kinetics  code  ANCON  (Ref.27).  For  the  purpose  of 
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Figure  3-7.  Benchmarking  of  DYNAM  Against  Analytical  Methods  - Positive 
Reactivity  Insertion  Case 


p(t)  = '0.00065 
p(0)  = 0.0 

N (0)  = 100.0  watts 
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Figure  3-8.  Benchmarking  of  DYNAM  Agaist  Analytical  Methods  - Negative 
Reactivity  Insertion  Case 
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comparison,  three  different  ramp  reactivity  insertion  rates  were 
selected.  Either  six  delayed  neutron  groups  or  one  delayed  neutron 
group  treatment  can  be  used  in  ANCON.  Figures  3.9  and  3.10  compare 
the  results  from  DYNAM  against  six  delayed  neutron  group  ANCON 
results  (note  that  the  one  delayed  group  ANCON  results  agree  exactly 
with  the  DYNAM  results). 

Since  ANCON  cannot  deal  with  circulating  fuel  type  reactors  (it 
assumes  that  the  fuel  is  stationary),  for  the  purpose  of  comparing  the 
results  of  DYNAM  with  those  of  ANCON,  the  fuel  circulation  effects  (for 
example,  the  loss  of  delayed  neutron  precursors  from  the  core,  the 
delayed  neutron  precursor  reintroduction,  etc.)  are  zeroed  in  DYNAM. 

Results 

The  developed  computer  program,  DYNAM,  is  now  used  to  study 
the  transition  to  burst  power.  As  mentioned  before,  the  program  is 
very  versatile  and  can  be  used  in  its  present  form  or  with  very  little 
modification  to  study  different  aspects  of  the  rapid  transition.  In  the 
following  sections,  a few  selected  capabilities  and  applications  of  the 
program  are  demonstrated. 

First,  two  separate  reactivity  insertions  are  applied  to  the  core 
to  bring  it  up  from  almost  zero  power  to  full  power  in  roughly  100 
seconds.  During  this  transition,  crucial  system  parameters  like  core 
pressure,  core  temperature,  etc.,  are  estimated  by  DYNAM.  Based  on 
the  behavior  of  these  parameters,  system  design  can  be  developed  or 
modified.  For  example,  if  one  finds  that,  for  some  desired  operating 
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Figure  3-9.  Benchmarking  of  DYNAM  Against  ANCON  - Case  1 
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Time  After  Ramp  Reactivity  Insertion  (sec) 


Figure  3-10.  Benchmarking  of  DYNAM  Against  ANCON  - Case  2 
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power  level,  the  average  core  temperature  is  higher  than  some 
acceptable  value,  the  system  design  can  be  modified,  e.g,,  to  increase 
the  fuel  flow  rate  through  the  core  so  as  to  reduce  the  average  core 
fuel  temperature.  This  increased  fuel  flow  rate  through  the  core  might 
require  changing  the  ratings  of  the  pumps  or  even  changing  the  core 
inlet  and  exit  nozzle  areas.  The  behavior  of  these  crucial  system 
parameters  can  also  be  used  in  selecting  the  best  possible  reactivity 
insertion  pattern  that  is  capable  of  taking  the  system  from  zero  power 
to  full  power  in  a safe  manner. 

The  way  DYNAM  is  set  up  currently,  to  bring  about  the  transition 
to  burst  power,  one  has  to  specify  the  desired  reactivity  variation  (this 
can  be  almost  any  conceivable  reactivity  variation).  The  program  Avill 
then  estimate  the  fuel  mass  variation  needed  within  the  core  to  bring 
about  this  reactivity  variation.  It  also  determines  the  fuel  flow  rates 
needed  in  order  to  obtain  the  necessary  fuel  mass  inside  the  core,  and 
the  core  inlet  pressure  is  automatically  adjusted  to  give  the  required 
fuel  flow  rates. 

There  are  many  types  of  reactivity  variations  that  can  be  applied 
to  take  the  core  from  zero  power  to  full  power.  For  the  purpose  of 
demonstration,  two  different  reactivity  insertions  are  selected.  The 
first  one  is  a pure  linear  (or  ramp)  reactivity  variation.  In  fact,  as 
shown  in  Fig.  3.11,  it  is  a combination  of  two  linear  reactivity 
variations.  To  simulate  the  conditions  in  an  actual  transition  scenario, 
the  core  is  initially  assumed  to  be  far  subcritical  (ptot  = -1.0).  This 


Figure  Not  to  Scale 


83 


..  8 


o 


o 

(N 


o 

o 

a 

o 

ed 


& 


tn 

o 

% 

o 


rt/w ) iC;iAnoB3H  Hcx)<ia 


FIG.  3.11  BPGCR  Input  Reactivity  Variation  With  Time  - Case  I 
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would  then  correspond  to  a situation  where  the  fuel  is  just  being 
introduced  into  a previously  unfuelled  (voided)  central  chamber. 

During  the  first  part  of  the  reactivity  variation— when  the  core  is 
subcritical,  a very  steep  reactivity  variation  of  0.02  Ak/k  per  second  is 
introduced  into  the  core  (for  a U235  fueled  system  with  a delayed 
neutron  fraction  of  about  0.0065,  this  number  translates  to  roughly 
$3.0  per  second).  Once  the  system  reactivity  becomes  zero  (keff  = 1.0), 

the  slope  of  the  linear  reactivity  variation  is  reduced  to  0.00015  Ak/k 
per  second  (roughly  $0.02  per  second). 

Figures  3.12  through  3.16  show  the  variation  of  some  of  the 
important  system  parameters  during  the  rapid  transition  to  about  100 
MWt  as  computed  by  DYNAM.  The  core  power  is  set  at  an  initial  value 
of  1000  watts  within  the  program.  This  might  correspond  to  a 
situation  where  the  bimodal  system  is  getting  ready  to  go  to  a high 
power  mode  of  operation  from  a low  power  mode,  and  some  power  is 
generated  within  the  BPGCR  due  to  neutrons  reaching  it  from  the 
surrounding  PGCRs.  The  behavior  of  the  core  power  (Fig.  3.12)  can  be 
explained  based  on  the  behavior  of  the  reactivity  variation  alone  as 
follows.  Initially,  when  the  core  is  far  subcritical,  the  core  power  falls 
as  time  increases  even  though  the  system  reactivity  is  increasing.  This 
happens  due  to  the  fact  that  the  system  is  so  far  subcritical  that  it  is 
not  able  to  maintain  a steady  state  neutron  population.  Then  the  total 
system  reactivity  reaches  a point  (still  subcritical)  where  the  rate  of 
fuel  mass  increase  within  the  core  forces  the  power  to  go  up  (note 
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that  the  core  power  is  given  by  an  expression  which  involves  a product 
of  the  flux  and  the  fuel  atom  densities  within  the  core.  So,  even 
though  the  flux  is  still  falling  (since  the  core  is  still  subcritical),  the 
rate  of  increase  of  the  fuel  atom  densities  compensates  for  this  effect 
and  the  net  result  is  a power  that  increases  with  time).  Then  at  50 
seconds  after  the  initiation  of  the  transition,  the  system  reactivity 
becomes  exactly  zero  (critical  stationary  system),  and,  at  this  point, 
within  DYNAM,  the  positive  reactivity  insertion  rate  is  reduced  to  a 
much  lower  value  (see  Fig.  3.11).  The  system  power  now  drops  slightly 
due  to  fact  that,  for  this  circulating  fuel  system,  there  is  some  loss  of 
delayed  neutrons  while  going  through  the  external  loop  and,  even 
though  the  stationary  system  reactivity  is  zero  at  time  t = 50  seconds, 
the  circulating  system  reactivity  is  still  negative  (system  multiplication 
factor  is  less  than  one).  Eventually,  at  t = 85  seconds,  the  circulating 
system  becomes  critical  and  further  reactivity  addition  leads  to  a rapid 
power  increase. 

To  understand  the  behavior  of  other  system  parameters  (like 
core  temperature,  core  pressure,  etc.)  the  effect  of  the 
balance -of-plant  also  has  to  be  taken  into  account.  The  effect  of 
reactivity  variation  alone  cannot  explain  the  behavior  of  these  system 
variables.  For  example.  Fig,  3.14  shows  that  the  core  temperature 
increases  from  the  very  beginning  of  the  transition  even  though  the 
core  power  is  dropping  during  the  initial  phases  of  the  transition.  The 
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Figure  3-12.  BPGCR  Thermal  Power  Vs.  Time  - Case  I 
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Figure  3-13.  BPGCR  Total  Fuel  (UF6+He)  Mass  Variation  With  Time  - Case  I 
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Figure  3-14.  BPGCR  Average  Core  Temperature  Variation  with  Time  - Case  I 
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Figure  3-15.  BPGCR  Average  Core  Pressure  Variation  with  Time  - Case  I 
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Figure  3-16.  BPGCR  Fuel  Mass  Outflow  Rate  Variation  with  Time  - Case  I 
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reason  for  this  behavior  of  core  temperature  can  be  directly  traced 
back  to  the  compressor  work  done  on  the  fuel/ working  fluid.  As 
mentioned  earlier,  since  this  is  a Brayton  cycle,  the  compressor  work 
needed  to  compress  the  gaseous  fuel  is  considerable.  This  effect 
dominates  during  the  initial  phase  of  the  transition  and  the  core 
power  (fission  power)  has  very  little  effect  on  the  core  temperature. 
The  work  done  by  the  compressor  depends  on  the  fuel  mass  flow  rate 
and,  hence,  initially,  during  the  rapid  insertion  of  reactivity  (addition 
of  fuel  to  the  core),  the  work  done  by  the  compressor  is  also 
increasing.  Only  towards  the  end  of  the  transition,  when  the  core 
power  becomes  significant  (above  a few  megawatts),  does  the  core 
power  starts  influencing  the  behavior  of  core  parameters  like  core 
pressure,  core  temperature,  etc..  From  Fig.  3.16,  we  can  see  that  most 
of  the  fuel  mass  flow  rate  increase  occurs  during  the  first  50  seconds 
(except  for  the  sharp  increase  towards  the  end)  of  the  transition.  This 
behavior  of  the  fuel  mass  flow  rate  gets  reflected  in  the  fuel 
temperature  profile— where  the  core  temperature  pretty  much  comes 
to  a steady  state  after  the  first  50  seconds  or  so.  The  discontinuity  in 
the  core  parameters  around  50  seconds  after  the  initiation  of  the 
transition  can  be  attributed  to  the  discontinuity  in  the  reactivity 
insertion  rate  at  this  point  in  time  (i.e.,  a change  in  the  reactivity 
insertion  rate  results  in  a change  in  the  fuel  mass  inflow  rate  and, 
hence,  a change  in  the  compressor  work  output).  The  core 
temperature  profile  does  not  show  any  sharp  discontinuity  at  50 


92 


seconds  because  of  the  damping  that  occurs  in  the  external  loop  (due 
to  the  heat  exchanger/radiator). 

The  type  of  reactivity  insertion  described  above  (pure  ramp 
insertion)  is  not  a good  one  to  employ  during  a real  transition.  As  can 
be  seen  from  Figures  3.12  through  3.16,  most  of  the  core  power 
increase  (and,  hence,  the  increases  in  the  core  temperature,  core 
pressure,  etc.)  occurs  very  rapidly  during  the  last  few  seconds  of  the 
transition.  This  means  that  some  type  of  drastic  controlling  action 
needs  to  be  taken  at  the  end  of  the  transition  period  to  prevent  power 
overshoot.  Also,  most  of  the  time,  the  safety  features  engineered  into  a 
reactor  are  designed  to  prevent  such  rapid  power  increases,  and  will 
cause  the  reactor  to  be  scrammed.  So.  instead  of  using  just  a pure 
ramp  insertion  of  reactivity,  a specific  reactivity  variation,  as  shown  in 
Fig.  3.17,  was  next  selected  for  the  purpose  of  demonstration.  This 
reactivity  variation  is  designed  to  produce  a fast  power  increase  to  a 
predetermined  constant  power  level.  As  shown  in  Fig.  3.17,  the 
reactivity  variation  is  a combination  of  two  linear  ramps  (with  different 
slopes)  and  an  exponentially  decaying  reactivity  variation.  By  adjusting 
the  peak  value  of  the  core  reactivity  (Pmax)  and  the  time  at  which  the 

exponentially  decaying  reactivity  variation  begins  (ti),  the  steady  state 

power  level  attained  can  be  varied  (see  Ref. 28  for  details  of  this  type  of 
reactivity  variation.  Ref. 28  also  has  details  about  other  types  of 
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Figure  3.17.  BPGCR  Input  Reactivity  Variation  With  Time  - Case  II 
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reactivity  variations  that  are  capable  of  producing  fast  monotonic 
power  increase  to  predetermined  constant  power  levels). 

Figures  3. 18  through  3.23  show  the  variation  of  the  most 
important  system  parameters  during  the  rapid  transition  as  predicted 
by  DYNAM.  Again,  the  behavior  of  the  system  parameters  can  be 
explained  as  a combined  effect  of  the  type  of  reactivity  variation 
introduced  into  the  core  during  the  transition,  plus  that  of  the 
balance-of-plant,  as  was  done  for  the  previous  case. 

Another  effect— which  has  some  significant  effect  on  the 
dynamics  of  the  system— can  be  noted  from  the  transition  scenario 
described  above.  The  beryllium  temperature  hardly  changes  during  the 
entire  transition.  This  is  the  case  even  when  there  is  no  heat  removal 
mechanism  provided  for  the  beryllium  moderator.  In  the  program,  the 
beryllium  temperature  is  initialized  as  800  K.  This  might  represent  a 
case  where  the  beryllium  temperature  is  maintained  at  800  K during 
the  low  power  mode  of  operation  (see  discussion  in  Chapter  1).  The 
sluggish  behavior  of  the  beryllium  moderator  temperature  is  due  to  the 
large  thermal  inertia  of  the  moderator.  One  outcome  of  this  sluggish 
behavior  of  the  beryllium  temperature  is  that  the  reactivity  feedback, 
due  to  the  moderator  temperature  change,  can  be  safely  neglected 
during  the  rapid  transition.  It  needs  to  be  taken  into  account  only  for 
dynamic  analysis  of  prolonged  operational  transients  at  high  power 
levels. 
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Figure  3-19.  BPGCR  Fuel  (UF6+He)  Mass  Variation  with  Time  - Case  II 
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Figure  3-20.  BPGCR  Average  Core  Pressure  Variation  with  Time  - Case  II 
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Figure  3-2 1 . BPGCR  Average  Core  Temperature  Variation  with  Time  - Case  II 


180 


99 


o 

V 

0) 


a 

o 


0< 

:3 

I 

■M 

(/) 

U 

4) 


4) 

B 

H 


(s/^H)  Momno  ssbi\i  lan^  ajoo 


Figure  3-22.  BPGCR  Fuel  (UF6+He)  Mass  Outflow  Rate  Vs.  Time  - Case  II 
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Figure  3-23.  BPGCR  Average  Moderator  Temperature  Vs.  Time  - Case  II 
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The  predictions  of  DYNAM  during  the  transition  can  be  very 
helpful  in  designing  the  system  or  modifying  the  design  of  the  system. 
An  example  of  this  capability  of  DYNAM  can  be  seen  by  considering 
Figure  3.21.  Here  the  variation  of  the  average  core  temperature  is 
given  when  the  reactivity  insertion  shown  in  Fig.  3.17  is  applied  to  the 
system.  Now  suppose  it  is  decided  that  the  steady  state  (at  the  end  of 
the  rapid  transition)  mean  core  temperature  is  not  high  enough.  One 
way  of  increasing  this  average  temperature  (without  increasing  the 
core  power)  is  to  decrease  the  fuel  mass  flow  rate  through  the  core. 
This  effect  is  shown  in  Fig.  3.24— where  the  core  average  temperature 
variation  is  given  for  the  same  reactivity  insertion  shown  in  Fig.  3.17, 
but,  with  a core  exit  throat  area  of  0.03  m2.  For  the  cases  shown  in 
Figures  3.18  through  3.23,  the  throat  area  of  the  converging- diverging 
nozzle  at  the  core  exit  is  fixed  at  0.05  m2.  This  area,  to  some  extent, 
dictates  the  fuel  mass  flow  rate  through  the  core  (see  Eq.3.10).  (Note 
that  at  the  core  exit,  the  flow  is  kept  choked  so  as  to  obtain 
supersonic  velocities  at  the  MHD  duct  inlet). 

It  can  be  seen  that  the  steady  state  core  temperature  has  now 
increased  to  about  2400  K from  the  previous  value  of  about  2000  K. 
Figure  3.25  shows  the  new  fuel  mass  outflow  rate  from  the  core. 
Comparing  Fig.  3.25  with  Fig.  3,22,  it  can  be  seen  that  the  fuel  mass 
flow  rate  is  now  lower  (resulting  in  a higher  core  temperature). 
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Figure  3-24.  BPGCR  Average  Core  Temperature  Vs.  Time 
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Figure  3-25.  BPGCR  Fuel  (UF6+He)  Mass  Outflow  Rate  Vs.  Time 


CHAPTER  4 

BPGCR  STABILITY  ANALYSIS 
Introduction 

Investigation  of  the  inherent  stability  constitutes  an  essential 
part  of  a reactor  design  process.  This  is  particularly  true  for  a new 
type  of  reactor  system  such  as  the  BPGCR  concept  being  investigated 
in  this  study.  A fundamental  consideration  in  such  an  investigation  is 
to  determine  whether  the  system  possesses  inherent  self-destruction 
tendencies,  i.e.,  to  determine  whether  a given  djmamical  system  (in 
this  case,  the  reactor)  will  eventually  return  to  an  equilibrium  state 
after  a finite  perturbation.  A reactor  is  clearly  unstable  if  a reactivity 
perturbation  can  cause  an  unlimited  increase  in  neutron  flux  or  power 
output.  This  is  true  whether  the  increase  in  power  or  neutron  flux 
occurs  continuously  or  as  a series  of  oscillations  of  increasing 
amplitude.  For  the  system  to  be  considered  inherently  stable  (i.e., 
without  the  help  of  any  external  controlling  action),  a perturbation 
during  steady  state  operation  must  eventually  lead  to  the  same  or 
another  equilibrium  state  of  the  system.  The  transient  behavior  in 
going  from  one  steady  state  to  another  should  also  be  satisfactory,  i.e., 
the  power  output  (or  flux)  should  not  undergo  a series  of  oscillations  of 
unacceptable  magnitude  before  settling  down. 
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Several  approaches  may  be  used  for  stability  analysis  of  a reactor 
system.  A complete  study  of  power  reactor  dynamics  would  take  into 
account  the  inherent  nonlinearity  of  the  reactivity  feedback.  This 
nonlinearity  is  due  to  the  dependence  of  the  system  reactivity  on  core 
power.  The  transient  response  of  a nonlinear  system  can  be  obtained 
by  solving  simultaneously  the  system  equations.  As  an  example  of  this 
procedure,  using  the  computer  program,  DYNAM  (described  in 
Chapter  3),  the  behavior  of  the  BPGCR  after  a reactivity  perturbation 
can  be  studied.  DYNAM  numerically  solves  the  nonlinear  neutron 
kinetics,  thermodynamics,  and  heat  transfer  equations  of  the  system. 

The  power  output  versus  time,  as  predicted  by  DYNAM,  for  a 
single  chamber  BPGCR  system,  initially  operating  at  a steady  state 
power  of  roughly  100  MWt,  following  a one  dollar  ($1)  external 
positive  step  reactivity  insertion  is  shown  in  Fig.  4. 1 . Note  that  for  a 
circulating  fuel  reactor,  the  reactivity  which  corresponds  to  a dollar  is 
not  the  same  as  that  for  a stationary  fuel  reactor.  Since  the  effective 
delayed  neutron  fraction  is  reduced  in  a circulating  fuel  reactor  due  to 
the  loss  of  delayed  neutrons  that  are  emitted  external  to  the  core,  the 
reactivity  that  corresponds  to  a dollar  is  given  by  'a  J3'  , where  'a'  is 
called  the  delayed  neutron  attenuation  factor,  'a'  is  a function  of  both 
core  residence  time  and  external  loop  circulation  time  of  the  fuel  (see 
Eq.  4.4).  The  external  loop  circulation  time  is  held  constant  at  0.1 
second  for  this,  as  well  as  for  most  of  the  other  transients  considered 
in  this  chapter.  (In  an  actual  system,  the  external  loop  circulation  time 


Initial  Core  Power  = 100  MWt 
$1  Positive  External  Reactivity  Step 

(pext  = -t-0.0029) 
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Figure  4- 1 . BPGCR  Thermal  Power  Vs.  Time  After  a One  Dollar 
External  Positive  Reactivity  Insertion 
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depends  on  the  length  of  the  external  loop  and  the  average  fuel 
velocity  while  traversing  the  external  loop.  Typical  values  for  the 
external  loop  circulation  time  should  only  be  a few  tenths  of  a second, 
and,  for  realistic  values  of  external  loop  length  and  fuel  velocities,  a 
second  is  probably  an  upper  limit  for  the  external  loop  circulation 
time.)  The  core  fuel  residence  time  is  of  the  order  of  0.08  sec.  The 
system  reactivity  is  a combination  of  externally  applied  reactivity,  plus 
the  reactivity  due  to  the  fuel  mass  inside  the  core.  No  feedback  effects 
due  to  beryllium  moderator  temperature  change,  or  Doppler  effect  due 
to  fuel  temperature  change  are  considered  for  this  perturbation.  The 
beryllium  temperature  takes  a long  time  to  change  because  of  its  large 
thermal  inertia  and,  hence,  the  feedback  effect  due  to  this  change  is 
expected  to  be  negligible  for  transients  that  are  but  a few  seconds  in 
duration.  The  fuel  temperature  (Doppler  effect)  coefficient  of  reactivity 
is  neglected  because  of  its  small  magnitude  compared  to  the  fuel  mass 
(or  density)  feedback.  For  the  illustrated  transient,  the  inlet  stagnation 
pressure  is  kept  constant.  During  normal  operation,  the  core  inlet 
stagnation  pressure  is  constantly  monitored  and  is  adjusted  so  as  to 
obtain  the  necessary  fuel  mass  inside  the  core.  However,  for  this 
external  step  reactivity  insertion,  the  value  of  the  inlet  stagnation 
pressure  is  frozen  at  its  initial  steady  state  value,  and  is  kept  constant 
during  the  transient. 

As  shown  in  the  figure,  as  soon  as  the  external  step  reactivity  is 
applied,  the  core  power  begins  to  rise.  This  results  in  an  increase  in 
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the  core  temperature  and  pressure.  But  since  the  core  inlet  stagnation 
pressure  is  kept  constant,  the  increase  in  core  pressure  has  two 
effects;  first  the  fuel  mass  flow  rate  at  the  core  inlet  starts  decreasing 
(the  core  inlet  is  always  maintained  at  non-choked  flow  conditions) 
and  second,  at  the  same  time  the  fuel  mass  flow  rate  at  the  core  exit 
starts  increasing.  These  two  effects  drive  the  total  core  reactivity 
down  and  the  core  power  starts  to  decrease.  The  rate  of  power 
decrease  is  halted  (or  reduced)  only  when  the  core  pressure  has  fallen 
enough  to  make  the  constantly  held  inlet  stagnation  pressure  high 
enough  (with  respect  to  the  core  pressure)  to  make  the  total  core 
reactivity  (fuel  mass)  go  up  again. 

It  can  be  seen  from  the  figure  that  the  strong  fuel  mass  (or 
density)  feedback  has  the  effect  of  bringing  the  system  under  control 
rather  quickly  (in  less  than  a second)  even  after  a hefty  reactivity 
perturbation  ($1  step  insertion).  The  amplitude  of  oscillation  is  also 
very  modest  (less  than  50%  of  the  initial  steady  state  value)  during  the 
initial  phase  of  the  transient.  The  system  behaves  in  a safe  manner 
after  the  reactivity  perturbation,  i.e.,  no  runaway  power  output  ensues 
after  the  perturbation,  even  though  no  external  controlling  action  is 
provided  (except  for  maintaining  the  core  inlet  stagnation  pressure 
constant  during  the  transient).  The  BPGCR  internal  reactivity 
(reactivity  due  to  the  fuel  mass  inside  the  core)  is  plotted  versus  time 
for  this  transient  in  Fig.  4.2.  During  the  entire  transient,  the  externally 
applied  reactivity  is  held  constant  at  one  dollar. 


0.0040 


109 


( n/yiv)  iBtua^ui  ajoQ 


Figure  4-2.  BPGCR  Internal  Reactivity  Vs.  Time  After  a One  Dollar 
External  Positive  Reactivity  Insertion 
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From  Fig.  4.2,  it  can  be  seen  that  the  net  (internal  plus  external) 
system  reactivity  remains  positive  during  the  entire  2 second 
transient,  and  yet,  the  system  power  is  decreasing.  This  is  happening 
because  of  the  loss  of  the  delayed  neutrons  due  to  fuel  circulation,  i.e., 
the  decay  of  some  of  the  delayed  neutron  precursors  while  in  the 
external  loop  effectively  reduces  the  delayed  neutrons  available  to  the 
system  for  multiplication.  So,  the  circulating  fuel  reactor  reactivity  has 
to  be  positive  (and  not  just  zero)  for  the  system  power  to  remain 
steady  or  to  increase. 

From  Fig.  4.1,  it  can  be  seen  that  an  external  positive  reactivity 
step  results  in  a lower  steady  state  power  compared  to  the  steady  state 
power  level  before  the  external  reactivity  insertion.  It  might  then  be 
inferred  that  a negative  external  reactivity  step  would  result  in  a 
steady  state  power  level  that  is  higher  compared  to  the  power  level 
before  the  reactivity  insertion;  this  would  be  a very  undesirable 
response.  So,  to  study  the  effect  of  a negative  reactivity  insertion,  a 
negative  one  dollar  reactivity  step  is  applied  to  the  BPGCR.  The 
resulting  power  behavior  is  given  in  Fig.  4.3.  It  can  be  seen  that,  for 
this  case  also,  the  final  steady  state  power  level  is  lower  than  the 
power  level  before  the  perturbation.  The  internal  reactivity  variation  is 
shown  in  Fig.  4.4.  Comparing  this  figure  with  Fig.  4.2,  it  can  be  seen 
that,  during  the  initial  oscillations,  the  behavior  of  the  internal 
reactivity  after  a negative  external  step  reactivity  insertion  is  almost  a 
mirror  image  of  that  after  a positive  external  reactivity  step  insertion. 
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Figure  4-3.  BPGCR  Thermal  Power  Vs.  Time  After  a One  Dollar 
External  Negative  Reactivity  Insertion 
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Figure  4-4.  BPGCR  Internal  Reactivity  Vs.  Time  After  a One  Dollar 
External  Negative  Reactivity  Insertion 
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This  response  of  the  BPGCR  is  quite  different  from  that  of  a 
conventional,  solid  core  (stationary  fuel)  reactor,  i.e,,  the  BPGCR 
steady  state  power  cannot  be  increased  merely  by  an  external  positive 
reactivity  addition.  Although  the  behavior  exhibited  in  Figs.  4.1  and  4.2 
may  be  a very  desirable  one  during  accident  situations  (abrupt 
insertion  of  positive  or  negative  reactivity),  it  can  pose  some  serious 
operational  problems  during  normal  reactor  operation.  During  normal 
reactor  operations,  a power  increase  can  be  brought  about  with  the 
help  of  methods  similar  to  those  employed  during  the  BPGCR  startup, 
i.e.,  a combination  of  external  positive  reactivity  addition  (e.g.,  by 
means  of  control  rods  or  drums),  plus  a change  in  the  core  inlet 
pressure.  This  method  is  discussed  in  detail  in  the  results  section  of 
Chapter  3. 

To  see  whether  the  delayed  neutron  circulation  effect  is 
influencing  the  system  power  behavior  significantly  during  reactivity 
transients,  two  additional  transients  are  studied;  in  the  first  one,  the 
external  loop  residence  time  of  the  fuel  is  increased  from  0.1  to  0.3 
seconds,  and  the  same  positive  external  reactivity  step  of  magnitude 
0.0029  Ak/k  is  applied  to  the  system  operating  at  a steady  state  power 
of  100  MWt.  In  the  second  one,  the  external  loop  circulation  time  is 
increased  from  0.1  seconds  to  0.5  seconds  and  an  external  positive 
reactivity  step  of  magnitude  0.0029  Ak/k  is  applied  to  the  system  (the 
core  fuel  residence  time  is  kept  fixed  at  0.08  seconds  in  both  cases).  A 
change  in  the  external  loop  circulation  time  affects  the  BPGCR 
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dynamics  by  changing  the  fraction  of  delayed  neutron  precursors  that 
are  reintroduced  into  the  core  (without  deca5dng)  after  traversing  the 
external  loop.  Figures  4.5  and  4.6,  respectively,  show  the  system 
power  behavior  for  the  above  two  cases.  For  the  purpose  of 
comparison,  the  results  from  the  transients  shown  in  Figs.  4.1,  4.5, 
and  4.6  are  combined  and  are  shown  together  in  Fig.  4.7. 

From  Fig.  4.7,  it  can  be  seen  that  the  variations  in  the  external 
loop  circulation  time  have  very  little  effect  on  the  initial  oscillatory 
behavior  of  the  system  (the  difference  in  amplitudes  of  the  maximum 
power  attained  by  the  system,  immediately  following  the  external  step 
reactivity  insertion,  is  less  than  5%  for  the  three  cases  shown)  but,  the 
changes  in  the  external  loop  circulation  time  do  significantly  affect  the 
final  power  attained  by  the  system  at  the  end  of  the  transient  (there  is 
a difference  of  roughly  25%  in  the  final  power  attained  by  the  BPGCR 
among  the  three  cases  shown  in  Fig.  4.7).  The  initial  oscillatory 
behavior  is  primarily  dictated  by  the  strong,  prompt  fuel  mass 
feedback  effect,  while  the  final  power  is  dictated  by  the  much  slower 
delayed  neutron  feedback  effect  due  to  the  fuel  circulation.  The  final 
powers  are  higher  for  the  cases  with  longer  external  loop  circulation 
times  (Ti  = 0,5  second  and  Ti  = 0.3  second)  compared  to  the  case  with 

the  shorter  external  loop  circulation  time  (li  = 0.1  second)  because  of 

the  fact  that,  although  the  same  (magnitude)  external  reactivity  step  in 
absolute  units  is  applied  to  the  system  in  all  cases,  for  the  longer 
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Figure  4-5.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Varying  External  Loop  Circulation  Time  - Case  1 
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Figure  4-6.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Varying  External  Lx)op  Circulation  Time  - Case  II 
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Figure  4-7.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 
Insertion  - Effect  of  Varying  External  Loop  Circulation  Time 
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external  loop  circulation  time  cases  (i.e.,  with  Ti  = 0.5  second  and  li  = 

0.3  second),  the  external  reactivity  step  corresponds  to  more  than  a 
dollar  in  relative  units.  For  a circulating  fuel  reactor,  the  magnitude  of 
relative  reactivity  that  corresponds  to  a dollar  is  given  by  ' a^',  where 
'a'  is  defined  as  the  delayed  neutron  attenuation  factor.  The  factor  'a' 
depends  on  the  core  residence  time  as  well  as  on  the  external  loop 
circulation  time  of  the  fuel  (see  Eq.  4.4  for  an  expression  for  'a').  The 
magnitude  of  the  applied  external  reactivity  step  corresponds  to 
exactly  a dollar  for  the  case  with  Xi  = 0.1  second  (and  Tc  = 0.08 

second).  This  difference  in  reactivity  insertion,  measured  in  relative 
units,  is  why  the  amplitude  of  oscillations  for  the  case  with  Xj  = 0.3 

second  is  slightly  higher  compared  to  the  case  with  Xi  = 0. 1 second, 

and  the  amplitude  of  oscillations  for  the  case  with  Xi  = 0.5  second  is 

slightly  higher  compared  to  the  cases  with  Xi  = 0. 1 second  and  Xi  = 0.3 
second. 

To  help  determine  the  error  associated  with  performing  a 
dynamic  analysis  of  the  BPGCR  with  only  one  delayed  neutron  group 
(as  opposed  to  performing  the  analysis  with  all  six  delayed  neutron 
groups),  a twenty  cent  (20 <t)  positive  reactivity  step  is  applied  to  the 
system  and  the  power  behavior  of  the  system  is  studied  for  the  cases 
with  one  delayed  neutron  group  and  with  six  delayed  neutron  groups. 
The  results  are  given  in  Figs.  4.8  and  4.9,  respectively.  It  can  be  seen 
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- One  Delayed  Neutron  Group  Case 
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Figure  4-9.  BPGCR  Thermal  Power  Vs.  Time  After  an  External 

Reactivity  Insertion  - Six  Delayed  Neutron  Group  Case 
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that,  for  the  examined  system,  neither  the  power  behavior 
immediately  folloAving  the  perturbation,  nor  the  final  steady  state 
power  attained  by  the  system  are  significantly  affected  by  inclusion  of 
all  six  delayed  neutron  groups. 

Comparing  Figs.  4.8  (one  delayed  neutron  group  case  for  20t 
positive  reactivity  step)  and  4.3  (one  delayed  neutron  group  case  for 
$1  negative  reactivity  step),  it  can  be  seen  that  the  final  steady  state 
power  attained  by  the  BPGCR  is  almost  the  same.  This  behavior— as 
mentioned  before— is  very  different  from  that  of  a conventional,  solid 
core,  stationary  fuel  reactor.  In  a conventional  reactor,  the  fuel  mass 
within  the  reactor  remains  constant  and,  if  a positive  (or  negative) 
external  reactivity  step  is  applied  to  such  a system  operating  at  steady 
state,  the  net  system  reactivity  will  remain  positive  (or  negative),  and 
the  system  power  will  correspondingly  increase  (or  decrease)  after 
the  step  reactivity  insertion.  For  the  BPGCR,  the  fuel  mass  within  the 
core  does  not  remain  constant  during  the  transient.  Immediately 
following  the  reactivity  perturbation,  the  amount  of  fuel  mass  within 
the  core  undergoes  oscillations  before  quickly  coming  to  an 
equilibrium.  The  magnitude  of  the  equilibrium  internal  reactivity 
(associated  with  the  fuel  mass  within  the  core)  is  such  that  the  total 
reactivity  (internal  plus  the  externally  applied)  will  bring  the  system  to 
a new  steady  state  power  level.  Note  that  because  of  the  loss  of  delayed 
neutron  precursors  from  fuel  circulation,  the  magnitude  of  the 
reactivity  needed  for  steady  state  power  is  greater  than  zero. 
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Since  the  fuel  mass  feedback  effect  is  dictating  the  system 
power  behavior  during  the  transients  described  above,  it  is  decided  to 
see  the  effect  of  changing  the  magnitude  of  the  fuel  mass  feedback 
coefficient  of  reactivity  on  the  system  behavior.  Under  normal 
conditions,  the  value  of  this  coefficient  changes  with  time  depending 
on  the  operating  conditions  (i.e.,  the  fuel  mass  within  the  core 
dictates  the  magnitude  of  this  coefficient).  To  see  the  effect  of 
changing  this  coefficient,  within  DYNAM,  the  magnitude  of  this 
coefficient  (Ap/AMf)  is  artificially  changed,  and  is  fixed  at  two  separate 

values  of  0.001  and  0.0005  (Ak/k  per  Kg  of  fuel),  respectively  (during 
normal  operating  conditions,  the  magnitude  of  this  coefficient  is  of  the 
order  of  0,035).  The  results  obtained  are  given  in  Figs.  4.10  and  4.11, 
respectively.  As  the  value  of  the  fuel  mass  feddback  coefficient  is 
reduced,  the  oscillations  previously  exhibited  by  the  system  (Figs.  4.1 
through  4.6),  disappear  completely.  This  is  an  expected  result  since, 
we  are  essentially  decreasing  the  negative  feedback  mechanism 
available  to  the  system.  The  fuel  mass  is  no  longer  as  strong  a feedback 
mechanism  as  it  was  before,  and  the  variations  in  core  pressure  and 
core  temperature  do  not  affect  the  system  reactivity  as  much  as  they 
did  previously.  But  if  we  reduce  the  fuel  mass  feedback  coefficient  too 
much,  the  system  becomes  unstable  and  the  power  starts  increasing 
without  limit  after  reactivity  perturbations  (see  Fig.  4,11). 
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Figure  4-10.  BPGCR  Termal  Power  Vs.  Time  After  an  External  Reactivity  Insertion  - 
Effect  of  Varying  the  Fuel  Mass  Feedback  Coefficient  - Case  I 
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Figure  4-11,  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity  Insertion  - 
Effect  of  Varying  the  Fuel  Mass  Feedback  Coefficient  - Case  II 
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At  the  other  extreme,  too  large  a value  for  this  coefficient  will 
also  make  the  system  unstable.  For  example,  results  for  the  case  where 
the  value  of  this  coefficient  is  increased  from  its  typical  value  of  0.035 
to  a value  of  0.1,  are  shown  in  Fig.  4.12.  The  oscillations  now  take 
more  time  to  die  down  compared  to  previous  cases,  and  it  is  found 
that  the  BPGCR  will  break  into  undamped  oscillations  if  the  value  of 
this  coefficient  is  increased  beyond  0.2.  An  optimum  value  of  this 
coefficient  can  be  selected  with  the  help  of  figures  like  the  ones 
shown  and  the  system  design  can  be  modified  so  as  to  operate  in  a 
regime  of  the  fuel  mass  versus  reactivity  curve  so  as  to  obtain  the 
desired  value  of  afm- 

The  coefficient  afm  is  obtained  from  the  slope  of  the  keff  versus 

core  fuel  mass  curve.  The  value  of  ttfm  can,  thus,  be  changed  by 

operating  at  a different  fuel  mass  loading.  In  order  to  maintain  the 
needed  keff,  this  will  require  some  form  of  reactivity  compensation. 
This  can  be  achieved  with  the  help  of  external  control  rods  or  drums 
located  in  the  external  moderator/reflector  region  or  by  adjusting  the 
gas  fuel  enrichment  or  by  adding  a burnable  poison  to  the  fuel  gas 
mixture. 

Finally,  to  the  see  the  effect  of  changing  the  inlet  pressure  to  the 
core,  at  the  same  time  that  the  one  dollar  reactivity  step  insertion  is 
applied  to  the  core,  the  value  of  the  core  inlet  stagnation  pressure  is 
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Figure  4-12.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity  Insertion  - 
Effect  of  Varying  the  Fuel  Mass  Feedback  Coefficient  - Case  III 
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raised  to  a much  higher  value  of  roughly  38  atm.  from  its  typical  value 
of  about  33  atm.  This  core  inlet  stagnation  pressure  is  kept  at  a 
constant  value  for  this  as  well  as  for  all  transients  described  previously. 
Since  the  fuel  mass  feedback  mechanism  is  dictating  the  system 
power  behavior  during  the  transient,  we  expect  any  effects  that  affect 
the  fuel  mass  inside  the  core  to  affect  the  power  behavior  also.  This 
can  be  seen  from  Fig.  4.13,  where  the  total  core  power  is  plotted  as  a 
function  of  time  after  the  initiation  of  the  transient.  Since,  now  we 
have  a higher  core  inlet  pressure,  more  fuel  is  forced  into  the  core 
and,  hence,  we  get  higher  power  oscillations  (due  to  higher  internal 
reactivity  oscillations).  Also,  the  frequency  of  oscillations  have 
increased  now  compared  to  the  previous  cases.  This  is  an  expected 
result  since  the  frequency  of  these  oscillations  depends  on  the  fuel 
transit  time  through  the  core  which  depends  on  the  fuel  inflow  and 
outflow  rates  which  are  in  turn  affected  by  the  higher  core  inlet 
stagnation  pressure. 

The  relative  magnitude  of  the  core  inlet  pressure  with  respect 
to  the  average  core  pressure  affects  the  core  fuel  loading  and,  hence, 
the  final  steady  state  power  level  attained  by  the  system  after  an 
external  reactivity  perturbation.  The  magnitude  and  sign  of  the 
externally  applied  reactivity  influences  only  the  initial  power 
oscillations  that  take  place  immediately  following  this  perturbation. 

The  effect  of  varying  the  core  inlet  stagnation  pressure  can  be  seen 
from  Figs.  4.14  through  4.16,  where  the  BPGCR  power  behavior 


Initial  Core  Power  =100  MWt 
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Figure  4-13.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Increasing  the  Core  Inlet  Stagnation  Pressure 
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Insertion  = 33  atm 

Constant  Core  Inlet  Pressure  = 34.1  atm 
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Figure  4-14.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Varying  the  Core  Inlet  Stagnation  Pressure  - Case  I 
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Figure  4-15.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Varying  the  Core  Inlet  Stagnation  Pressure  - Case  II 
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Figure  4-16.  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 

Insertion  - Effect  of  Varying  the  Core  Inlet  Stagnation  Pressure  - Case  III 


132 


following  a one  dollar  positive  reactivity  step  is  shown  for  three 
different  values  of  constant  core  inlet  stagnation  pressure.  It  can  be 
seen  from  the  figures  that,  as  the  value  of  the  core  inlet  stagnation 
pressure  is  increased  (relative  to  the  average  core  pressure),  the  final 
core  power,  as  well  as  the  amplitude  of  the  initial  power  oscillations, 
also  increase. 

For  all  the  reactivity  perturbations  examined  so  far,  the  core 
inlet  stagnation  pressure  is  kept  constant  during  the  transient.  It  is 
established  that  the  BPGCR  is  an  extremely  stable  configuration  under 
this  condition.  By  keeping  the  inlet  pressure  constant  during  the 
transient,  we  are,  in  effect,  providing  an  external  controlling 
mechanism.  To  see  whether  the  BPGCR  is  stable  even  without  this 
external  action  (i.e.,  to  see  whether  the  system  is  inherently  stable  or 
not),  a one  dollar  positive  external  reactivity  step  is  applied  to  the 
system  and  the  core  inlet  pressure  is  allowed  to  vary.  The  core  inlet 
pressure  at  any  given  time,  t,  is  set  equal  to  the  average  core  pressure 
which  had  existed  at  a time  Xi  back,  i.e.,  it  is  set  to  the  average  core 

pressure  at  time  (t-Ti),  where  Xi  is  the  external  loop  circulation  time. 

The  resulting  BPGCR  power  behavior  in  Fig.  4.17,  shows  that  the 
BPGCR  thermal  power  drops  more  than  two  orders  of  magnitude  in  a 
matter  of  few  seconds  (note  that  the  Y-axis  of  Fig.  4.17  has  a log  scale). 
So,  it  can  be  seen  that  when  the  core  inlet  pressure  is  allowed  to  float, 
the  BPGCR  becomes  even  more  stable  against  accidental  reactivity 
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Figure  4-17,  BPGCR  Thermal  Power  Vs.  Time  After  an  External  Reactivity 
Insertion  - Variable  Core  Inlet  Stagnation  Pressure  Case 
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insertions.  For  normal  reactor  operation,  this  power  behavior  is  much 
less  desirable  than  for  the  constant  core  inlet  pressure  case  because  of 
the  large  drop  in  power  Avithin  a short  period  of  time,  as  well  as  the 
fact  that  the  power  level  continues  to  drop  without  coming  to  a new 
steady  state  power  level  (not  shown  in  Fig.  4.17). 

The  above  described  method  of  analyzing  the  transient  response 
(stability)  of  a reactor  by  solving  simultaneously  the  set  of  nonlinear 
system  equations  can  rapidly  become  impractical  since  the  complete 
reactor  system  is  quite  complex,  and  usually  many  equations  are 
required  to  accurately  describe  the  behavior  of  the  system,  and  many 
transients  for  many  different  values  of  system  parameters  need  to  be 
examined  to  map  the  entire  operating  range.  Solving  these  nonlinear 
system  equations  can  be  both  time  consuming  as  well  as  tedious. 

A second,  and  less  complicated,  approach  to  studying  the 
stability  is  to  linearize  the  feedback  terms  in  the  system  equations  and 
use  the  well-developed  methods  of  linear  feedback  theory  for  stability 
analysis.  Some  important  conclusions  about  a nonlinear  system  can 
usually  be  made  by  stud5ang  its  associated  linear  approximation.  In 
particular,  the  stability  of  a nonlinear  system  for  small  perturbations 
may  often  be  deduced  by  examining  the  stability  of  an  associated  linear 
system. 

This  type  of  analysis  is  particularly  useful  for  obtaining  a good 
understanding  of  the  general  nature  of  the  reactor  behavior.  Effects  of 
changes  in  input  reactor  parameters  can  be  readily  identified,  and 
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regions  of  stability  can  be  mapped  for  various  ranges  of  reactor 
parameters.  Such  a linear  stability  analysis  of  the  BPGCR  is  undertaken 
in  the  following  section.  The  methods  of  analysis  include  the  familiar 
frequency  response  methods  (Bode  plots,  Nyquist  plots,  etc.)  and  pole 
configuration  (root  locus)  methods.  Details  of  these  methods  are  not 
given  in  this  dissertation,  but  can  be  found  in  a variety  of  system 
dynamics  textbooks  including  Ref,  26. 

BPGCR  Linear  Stability  Analysis 

Dynamically  the  most  important  characteristics  of  the  gas  core 
reactor  being  investigated  in  this  study  are  that  the  fuel  circulates  and 
that  the  fuel  is  highly  compressible.  The  fuel  circulation  acts  to  reduce 
the  effective  delayed  neutron  fraction  (in  general,  the  lower  the  value 
of  delayed  neutron  fraction,  the  more  difficult  it  becomes  to  control 
the  reactor— a reason  why  plutonium  fueled  reactors  are  more  difficult 
to  control  than  uranium  fueled  ones).  From  a control  standpoint,  the 
delayed  neutron  fuel  circulation  effect  acts  as  a positive  feedback  (i.e., 
if  the  reactor  power  goes  up,  the  delayed  neutrons  reintroduced  into 
the  core  after  going  through  the  external  loop  will  also  go  up,  and  this 
tends  to  drive  the  reactor  power  up  further)  while  the  gas  density 
effect  is  a negative  feedback.  These  are  two  opposite  (in  sign)  but 
powerful  feedback  control  mechanisms.  Fuel  circulation  also  acts  to 
reduce  the  rate  of  fuel  temperature  change  during  a power  transient. 

For  a linearized  system,  the  transfer  function  is  defined  as  the 
Laplace  transform  of  a selected  output  variable  divided  by  the  Laplace 
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transform  of  a selected  input  variable.  Its  value  is  a complex  number 
describing  both  amplitude  and  phase  relationships.  For  a reactivity 
input,  the  reactor  transfer  function  is  the  ratio  of  the  Laplace 
transform  of  the  steady  state  flux  response  to  the  Laplace  transform  of 
the  reactivity  input,  i.e.,  for  an  input  reactivity  p(t)  = po  Sin  (Ot,  the 

flux  will,  eventually  (after  the  transients  have  died  down),  be  given  by 

<t>  (t)  = 00  sin  (cot  + 0) 
and  the  transfer  function  is  given  by 

G Oco)  = ^ exp  (j0) 

PO 

where  '0o'  is  the  amplitude  of  the  of  the  flux  response,  'po'  is  the 

amplitude  of  the  reactivity  input,  and  '0'  is  the  phase  shift. 

The  open  loop  transfer  function  is  the  response  of  the  reactor  in 
the  absence  of  any  feedback  (temperature  coefficients  of  reactivity, 
etc.).  This  transfer  function  is  characteristic  of  a reactor  operating  at 
zero  power  where  the  flux  level  is  so  low  that  changes  in  the  flux  level 
do  not  significantly  affect  the  temperature  or  other  feedback  variables. 
The  closed  loop  transfer  function  is  the  response  of  the  reactor  at 
power,  including  all  feedback  mechanisms,  whereby  the  power  level 
affects  the  reactivity  of  the  reactor.  References  26  and  29  contain 
more  details  about  transfer  functions,  and  how  they  can  be  derived 
from  the  system  equations,  and  how  they  can  be  employed  to  obtain 
information  about  the  system  stability. 
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If  we  consider  an  elementary  critical  reactor  (i.e.,  one  with  no 
feedback  and  no  fuel  circulation),  the  transfer  function  of  this  reactor 
can  be  easily  obtained  from  the  neutron  kinetics  equations  (see  Ref. 
26),  and  is  given  below  for  the  case  of  one  delayed  neutron  group. 

G(s)  = ^ (s  + A.) 

A s (s+ p/A)  (4-1) 

where 

no  = equilibrium  value  of  the  neutron  density  (#/cc), 

s = Laplace  transform  variable, 

X = average  value  of  decay  constant  for  the  delayed  neutron 
precursors  (1/sec), 

A = neutron  generation  time  (sec),  and 

3 = total  delayed  neutron  fraction. 

If  we  analyze  this  transfer  function  of  an  elementary  reactor 
using  any  of  the  standard  linear  stability  analysis  tools,  we  will  find  that 
this  reactor  is  unstable.  For  example.  Fig.  4.18  shows  the  Bode  plot  of 
this  transfer  function.  A Bode  plot  is  a simple  plot  of  the  magnitude 
and  phase  angle  of  the  transfer  function  versus  frequency  (often.  Bode 
plots  are  used  to  display  experimental  results  in  a form  from  which  an 
approximate  transfer  function  can  readily  be  derived).  Figure  4.18  (in 
which  only  the  amplitude  plot  is  given  for  the  sake  of  illustration) 
shows  that  the  reactor  becomes  unstable  at  zero  frequency,  i.e,,  the 
reactor  has  infinite  gain  at  zero  frequency.  Physically,  this  means  that 


l.OOE+13 


138 


Figure  4-18.  Bode  Plot  for  Reactivity  Transfer  Function  of  an  Elementary, 
Zero  Power,  Non-Circulating  Fuel  Reactor 
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when  a small  amount  of  positive  reactivity  is  inserted  into  a critical 
reactor,  its  power  level  rises  indefinitely.  For  this  Bode  plot,  'no'  is 

taken  to  be  equal  to  1.0  E+8,  and  the  rest  of  the  parameters  have  the 
same  value  as  were  used  in  DYNAM  (i.e.,  a delayed  neutron  fraction  of 
0.0065,  a neutron  generation  time  of  6.0E-4  sec.,  and  an  average 
precursor  decay  constant  of  0.08  sec-i). 

In  a real  reactor,  whenever  operating  at  power,  there  are  always 
feedback  effects  present  which  tend  to  alter  the  response  of  the 
reactor  to  external  inputs  or  perturbations.  When  feedback  effects  are 
present,  the  elementary  reactor  transfer  function  is  modified  as  shown 
in  Fig.  4.19,  and  the  overall  reactor  transfer  function  is  given  by 

Gtot(s)  = G(s)  / (1  + G(s)H(s))  (4-2) 

where 

G(s)  = elementary  zero  power  reactor  transfer  function,  and 

H(s)  = feedback  transfer  function 

If  we  now  include  the  effect  of  fuel  circulation  on  the  elementary 
zero  power  reactor  model,  we  have  to  modify  the  elementary  reactor 
transfer  function  given  above.  The  circulating  delayed  neutron 
precursors  tend  to  act  like  a delayed  positive  feedback  on  the  reactor 
transfer  function.  To  understand  this  phenomenon,  consider  a 
circulating  fuel  reactor  operating  at  steady  state.  Of  all  the  delayed 
neutrons  precursors  produced  in  the  core,  only  a certain  fraction  will 
decay  within  the  core  before  leaving  the  core.  Some  of  the  precursors 
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will  decay  while  they  are  still  in  the  external  loop  and,  hence,  will  not 
contribute  to  the  core  neutron  population.  The  ones  that  survive 
without  decaying  while  traversing  the  external  loop  will  reenter  the 
core  and  then  can  decay.  When  the  reactor  is  operating  at  steady  state, 
a dynamic  equilibrium  is  established  for  the  delayed  neutron 
population  within  the  core.  Now,  if  for  some  reason,  the  core  power 
goes  up,  the  delayed  neutron  precursor  density  will  also  increase.  This 
increase  in  the  precursor  density  manifests  itself  as  an  increase  in  the 
fraction  of  precursors  that  are  reintroduced  into  the  core  after  some 
delay  (depending  on  the  external  loop  circulation  time)  and,  hence, 
will  tend  to  drive  the  power  further  up.  The  reverse  of  this  process 
occurs  when  the  power  goes  down. 

Using  the  equations  given  in  Ch.3  (Eq.  3.12)  for  the  precursor 
density  in  a circulating  fuel  reactor,  one  can  derive  the  transfer 
function  for  a zero  power,  circulating  fuel  reactor  (see  Ref.  29  for 
details  of  the  derivation).  It  is  given  by  the  following  expression: 

G(s)  = M£)  = 

s2I  + s ja  + + apj  + XP(a-l)  + 

where 

C2  = 1 - e-^'^i  e ®^i 

Xc  = mean  core  fuel  residence  time  (sec). 
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Xi  = mean  external  loop  eirculation  time  (sec), 

X = average  decay  constant  for  delayed  neutron  precursors 
(sec-i), 

= total  delayed  neutron  fraction,  and 
a = delayed  neutron  attenuation  factor  (see  discussion  below) 

In  the  above  equation,  the  quantity  "a"  is  called  the  delayed 
neutron  attenuation  factor.  Its  effect  is  to  reduce  the  delayed  neutron 
fraction  in  a circulating  fuel  reactor  by  a factor  "a".  We  can  intuitively 
infer  that  the  effective  delayed  neutron  fraction  in  a circulating  fuel 
reactor  has  to  be  less  than  that  for  a non-circulating  fuel  case.  Also, 
intuitively,  we  can  say  that  this  factor  "a"  has  to  depend  on  the  relative 
times  the  precursors  spend  within  the  core  and  in  the  external  loop. 
For  example,  if  the  external  loop  circulation  time  is  large,  a greater 
fraction  of  delayed  neutron  precursors  will  decay  while  in  the  external 
loop.  We  can  formally  derive  the  expression  for  "a",  to  obtain 

a = ^ 

Xxc  + i-e-^^i  (4-4) 

(See  Ref.  30.  for  details  of  this  derivation.) 

For  the  gas  core  reactor  being  investigated,  in  the  program 
DYNAM,  the  external  loop  circulation  time  is  generally  fixed  at  0. 1 
second.  The  core  residence  time  of  the  fuel  depends  on  the  velocity  of 
fuel  within  the  core  but  is  generally  in  the  range  of  a few  tens  of 
milliseconds  (20-40)  corresponding  to  fuel  velocities  of  the  order  of 
50  to  100  m/sec  (assuming  a core  length  of  about  2 meters).  The  Bode 
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plot  (magnitude)  for  this  zero  power,  circulating  fuel  reactor  transfer 
function  is  given  in  Fig,  4.20  (a  5th  order  Fade  approximation  is  used 
to  represent  the  delay  terms  of  Equation  (4,3)).  Comparing  this  figure 
with  Fig. 4. 18  (which  shows  the  transfer  function  for  a non-circulating 
fuel  reactor),  it  can  be  seen  that  the  gain  is  higher  (i.e.,  the  system  is 
less  stable)  at  lower  frequencies  for  the  circulating  fuel  reactor— an 
expected  result  due  to  the  lower  effective  value  of  the  delayed  neutron 
fraction  for  the  circulating  fuel  reactor. 

Inherent  Stability  of  the  BPGCR  Attached  to  an  External  Loop 

In  the  previous  section,  the  reactivity  transfer  function  of  the 
BPGCR  operating  at  zero  power  was  given.  However,  in  a reactor 
operating  at  power,  feedback  effects  are  always  present,  and  an 
external  loop  is  always  attached  to  the  reactor.  The  ultimate  concern 
is  the  stability  of  the  complete  system,  i.e.,  we  must  include  the  effects 
from  within  the  reactor,  as  well  as  the  effects  of  the  balance-of-plant 
(B.O.P),  to  check  for  the  overall  stability.  These  internal  feedbacks  and 
the  B.O.P  can  significantly  alter  the  response  of  a zero  power  reactor 
model  given  in  the  previous  section. 

From  an  elementary  transfer  function  standpoint,  the  complete 
system  can  be  considered  as  shown  in  Fig.  4.21.  Here  all  the  possible 
feedback  effects  that  might  affect  the  stability  of  the  system  are 
represented  as  two  distinct  loops.  Each  one  of  these  two  loops  can  be 
further  broken  down  to  include  component  transfer  functions. 
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Figure  4-20.  Bode  Plot  for  Reactivity  Transfer  Function  of  a Zero  Power, 
Circulating  Fuel  Reactor 
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Figure  4-2 1 . Simplified  Schematic  Diagram  of  the 
Complete  BPGCR  System 
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The  internal  loop  includes  the  various  power  coefficients  of 
reactivity,  i.e.,  functions  which  describe  how  the  reactor  power 
influences  the  reactivity  effects  due  to  changes  in  the  fuel 
temperature,  the  moderator  temperature,  the  fuel  mass,  etc..  The 
external  loop  represents  the  effect  of  components  like  the  heat 
exchanger,  the  radiator,  the  compressor,  etc.,  on  the  performance  of 
the  BPGCR.  These  effects  include  the  heat  transfer  from  the 
fuel/ coolant  to  the  heat  exchanger  (attenuation  of  power),  the  time 
involved  in  the  process  of  circulating  the  fuel  through  the  external 
loop  (time  delay),  etc.  The  combination  of  these  external  loop  transfer 
functions,  in  conjunction  with  the  zero  power  reactor  transfer 
function  and  its  associated  internal  feedback  functions,  determine  the 
stability  of  the  overall  system. 

External  Loop 

As  mentioned  before,  the  external  loop  is  associated  with  the 
B.O.P.  For  the  BPGCR,  the  fuel  and  the  working  fluid  (coolant)  are  one 
and  the  same.  The  gaseous  fuel/ working  fluid  mixture  leaving  the 
reactor  will  remove  most  of  the  heat  generated  within  the  core.  Under 
steady  state  conditions,  this  heat  is  given  up  by  the  fuel  gas  mixture 
while  going  through  the  external  loop  and  the  fuel  reenters  the  core  at 
a lower  temperature  after  some  delay.  The  amount  of  attenuation  (in 
temperature)  and  delay  the  fuel  gas  mixture  suffers  while  going 
through  the  external  loop  depends  on  the  components  of  the  external 
loop.  To  see  how  the  external  loop  can  affect  the  performance  of  the 
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reactor,  consider  a sudden  increase  in  the  power  generation  within 
the  core.  This  results  in  an  increase  in  the  temperature  of  the  coolant 
(fuel)  exiting  from  the  core  and,  hence,  that  of  the  fuel  gas  mixture 
reentering  the  core.  An  increase  in  these  two  temperatures  will  affect 
the  average  fuel  temperature  and  fuel  gas  mixture  density  within  the 
core  which,  in  turn,  can  influence  all  the  internal  feedback  loops  of 
the  reactor. 

To  derive  the  transfer  functions  for  the  external  loop,  first  the 
components  of  the  loop  have  to  be  selected  and  then  the  equations 
describing  their  performance  have  to  be  modelled.  Since  the  BPGCR  is 
only  in  a conceptual  stage  of  development,  there  is  a lot  of  flexibility  in 
selecting  the  various  components  of  the  B.O.P.  and  their  operating 
conditions  (for  example,  the  nature  of  the  working  fluid  in  the 
secondary  side  of  the  radiator/heat  exchanger  or  the  average 
secondary  fluid  temperature,  etc.).  For  the  purpose  of  this  analysis,  a 
simple  external  loop  as  shown  in  Fig.  4.22  is  selected.  It  consists  of  a 
turbine  (in  the  final  design  there  may  not  be  a turbine  as  such  but 
some  sort  of  an  energy  conversion  device  such  as  an  MHD  generator 
will  be  there— which  is  represented  here  as  a turbine),  a heat 
exchanger /radiator  combination  for  rejecting  the  waste  heat,  and 
finally  a compressor  for  pumping  the  fuel  back  to  the  core  at  the 
appropriate  pressure. 

The  equations  describing  each  of  the  above  components  are 
given  in  Appendix  B.  Since  this  is  a linear  stability  analysis,  only  small 
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variations  around  some  mean  operating  conditions  are  considered,  i.e., 
all  parameters  are  treated  as  perturbation  quantities  or  deviations 
from  their  steady  state  values.  The  transfer  function  representation  of 
the  external  loop  can  be  represented  as  shown  in  Fig.  4.23  by  using 
the  equations  given  in  Appendix  B.  As  shown  in  the  figure,  the  total 
time  delay  in  the  external  loop  is  handled  by  two  lumped  delays  (G6)— 
one  for  the  transfer  of  fuel  gas  mixture  from  the  core  to  the  heat 
exchanger  and  the  other  for  the  fuel  gas  mixture  transfer  from  the 
heat  exchanger  back  to  the  core.  In  Fig.  4.23,  the  transfer  function  G5, 
represents  the  combined  effect  of  the  turbine,  the  heat 
exchanger/radiator,  and  the  compressor.  Gl,  G2,  G3  and  G4  represent 
the  transfer  functions  that  describe  how  the  core  outlet  temperature 
(T3')  is  affected  by  the  core  inlet  temperature  (T2  ),  core  thermal 

power  (N).  core  inlet  fuel  mass  flow  rate  (m2),  and  the  core  exit  fuel 

mass  flow  rate  (m3),  respectively.  Appendix  B gives  the  contents  of  the 

transfer  function  boxes  shown  in  Fig.  4.23. 

Internal  Feedback  Loon 

The  structure  of  this  loop  for  the  BPGCR  is  significantly  different 
form  that  of  a conventional,  solid  core  reactor.  Again,  the  main  reason 
for  this  difference  is  the  gaseous  nature  of  the  fuel  plus  the  fact  that 
BPGCR  is  an  open  flow  system.  For  the  BPGCR— unlike  for  a 
conventional  solid  core  PWR,  the  most  significant  component  of  the 
reactor  internal  feedback  loop  is  the  fuel  mass  (or  density)  coefficient 
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Figure.  4-23.  Transfer  Function  Representation  of  the 
BPGCR  External  Loop 
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of  reactivity.  In  a conventional  solid  core  reactor,  the  fuel  mass  within 
the  reactor  remains  constant  whereas  in  the  BPGCR,  the  gaseous  fuel 
mass  is  a strong  function  of  the  operating  conditions  (core  average 
temperature,  pressure,  etc.)  and,  hence,  constitutes  a strong  feedback 
mechanism.  The  other  components  of  the  internal  feedback  loop  like 
the  fuel  temperature  (Doppler)  coefficient  of  reactivity  and  the 
moderator  temperature  coefficient  of  reactivity  are  of  much  less 
importance  compared  to  the  fuel  mass  coefficient  of  reactivity. 

The  equations  describing  the  internal  feedback  loop  are  given  in 
Appendix  B.  As  a first  step  in  the  analysis  of  the  internal  feedback 
loop,  it  can  be  represented  as  shown  in  Fig.  4.24.  This  figure  shows 
the  three  major  feedback  loops  that  can  influence  the  stability  of  the 
BPGCR.  The  core  inlet  and  outlet  temperatures  combine  to  provide  a 
core  average  fuel  gas  mixture  temperature.  Changes  in  this  average  fuel 
gas  mixture  temperature  directly  gives  rise  to  one  component  of  the 
internal  feedback  loop,  viz.,  the  fuel  temperature  (Doppler)  coefficient 
of  reactivity.  The  fuel  gas  mixture  temperature  also  indirectly 
influences  the  moderator  temperature  and  the  fuel  mass  (density) 
components  of  the  internal  reactivity  feedback  loop  as  shown  in  the 
figure. 

Formulations  for  the  fuel  temperature  feedback  and  the 
moderator  temperature  feedback  parts  of  the  internal  feedback  loop 
are  relatively  straightforward.  The  average  core  fuel  temperature  is 
obtained  as  a simple  average  of  the  core  inlet  and  outlet  temperatures. 
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Figure  4-24.  Transfer  Function  Representation  of  the  BPGCR  Internal 
Reactivity  Feedbacks 
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and  the  fuel  temperature  coefficient  of  reactivity  is  estimated  using 
123-group  XSDRNPM  (Sn  transport  theory)  calculations  (see  Chapter 

3.  Tables  3.5  and  3.6).  For  the  moderator  temperature  coefficient 
component,  the  equations  relating  the  variation  of  the  average 
moderator  temperature  as  a function  of  the  variation  of  the  average 
fuel  temperature  are  given  in  Appendix  B.  The  moderator  temperature 
coefficient  of  reactivity  is  also  estimated  using  123-group  XSDRNPM 
calculations  (see  Chapter  3,  Table  3.7). 

The  analysis  of  the  gaseous  fuel  mass  (or  density)  component  of 
the  internal  feedback  loop  is  more  involved  and  has  to  be  performed 
in  steps.  The  fuel  mass  inside  the  core  is  affected  by  a number  of 
variables  including  the  core  pressure,  core  temperature,  core  inlet 
stagnation  pressure,  fuel  mass  inflow  and  outflow,  etc..  The  complete 
set  of  equations  interrelating  all  these  variables  are  given  in  Appendix 
B.  Rearrangement  of  these  equations  and  some  algebraic  manipulations 
(see  Appendix  B),  yield  the  block  diagram  given  in  Fig.  4.25--which 
shows  the  relationship  among  the  involved  variables. 

For  the  initial  stability  analysis,  we  will  consider  only  the  fuel 
mass  component  of  the  internal  feedback  loop.  This  is  the  strongest  of 
the  three  internal  feedback  paths.  The  fuel  mass  feedback  effect  can 
be  added  to  Fig.  4.25  since  it  gives  the  variation  of  fuel  mass  within  the 
BPGCR.  Adding  the  fuel  mass  feedback  coefficient  to  this  fuel  mass 
variation,  yields  the  variation  of  that  component  of  internal  reactivity 
which  is  due  to  the  fuel  mass  variation.  Adding  a zero  power 
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Figure  4-25.  Transfer  Function  Representation  of  the  BPGCR  Fuel  Mass  Feedback 
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circulating  fuel  reactivity  transfer  function  to  this  internal  reactivity 
variation  gives  the  thermal  power  variation  within  the  BPGCR.  This  is 
shown  in  Fig.  4.26— where,  in  addition  to  the  fuel  mass  feedback 
component,  the  transfer  functions  which  represent  the  influence  of 
core  inlet  temperature,  core  thermal  power,  core  inlet  and  exit  fuel 
mass  flow  rates  on  the  core  exit  temperature  (see  Fig.  4.23),  are  also 
included.  Figure  4.26  now  represents  the  complete  BPGCR  except  for 
the  balance-of-plant.  Before  the  stability  analysis  of  the  system  can  be 
performed,  the  mass  feedback  branch  of  the  internal  feedback  loop 
given  above  has  to  be  further  reduced.  After  reducing  the  block 
diagram  given  in  Fig.  4.26,  we  can  combine  it  with  the  balance-of-plant 
loop  given  earlier  (Fig.  4.23)  to  obtain  a simplifled  representation  of 
the  complete  BPGCR  system  as  shown  in  Fig.  4.27.  This  is  the  model 
that  is  used  for  the  stability  analysis  of  the  system. 

The  linear  stability  analysis  of  the  BPGCR  is  performed  using  a 
computer  program  called  CC.  It  is  a personal  computer  based 
commercial  software  package  which  is  capable  of  performing  both 
classical  control  theory  analysis,  as  well  as  state  variable  analysis. 

For  performing  a classical  stability  analysis  (frequency  response, 
root  locus,  etc.)  using  CC,  the  complete  system  given  in  Fig.  4.27  has 
to  be  reduced  to  a single  transfer  function,  i.e.,  a single  block 
representing  the  relationship  between  an  input  and  an  output  as 


shown  below: 
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Figure  4-26.  Transfer  Function  Representation  of  the  Complete  BPGCR  with  Fuel 
Mass  Feedback  Excluding  the  Balance-of-Plant 
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Figure  4-27.  Simplified  Transfer  Function  Representation  of  the  Complete  BPGCR 
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Input 


Output 


To  get  the  system  block  diagram  into  such  a form,  first  the  external 
loop  can  be  eliminated.  To  do  this,  we  have  to  derive  two  transfer 
functions:  one  representing  the  variation  of  the  core  outlet 
temperature,  T3',  with  respect  to  the  variations  in  the  core  power 

(3T3  / 3N),  and  a second  one  representing  the  variations  of  the  core 
inlet  temperature,  T2',  with  respect  to  the  variations  in  the  core  power 
OT2  /3N).  These  two  transfer  functions  can  be  derived  easily  from  the 
external  loop  representation  given  in  Fig.  4.27  (the  derivation  is  given 
in  Appendix  B).  With  the  help  of  these  two  transfer  functions,  the 
block  diagram  of  Fig.  4.27  can  be  reduced  to  the  one  shown  in  Fig. 

4.28.  Fig.  4.29  shows  the  same  overall  transfer  function  block  diagram 
in  a more  conventional  manner,  where  an  external  reactivity  input  is 
also  shown.  This  external  reactivity  could,  for  example,  represent  the 
control  rods  or  an  adjustable  external  reflector,  or  some  type  of  an 
accidental  reactivity  insertion.  The  total  system  reactivity  is  given  by 
the  sum  of  this  external  reactivity  and  the  internal  reactivity. 

The  block  shown  in  Fig.  4.29  can  now  be  used  for  obtaining  the 
overall  system  transfer  function,  9N  / 3pext.  as 
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Figure  4-28.  Reduced  Transfer  Function  Representation  of  the 
Complete  BPGCR 


Figure  4-29.  Reduced  Transfer  Function  Representation  of  the 
Complete  BPGCR  in  Standard  Form 


dN  ^ ZPR 

dpext  1-  ZPR  am  (G22  G24  + G21  G23) 
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(4-5) 

where  ZPR  is  the  zero  power,  circulating  fuel  reactor  transfer 
function. 

Linear  Stability  Analysis  of  the  Overall  System  Transfer  Function 

Equation  4,5  is  of  the  familiar  closed  loop  transfer  function  form 
Y(s)  = G(s)  / (1  + G(s)H(s)) 

where  G(s)  represents  the  feed  forward  transfer  function  and  H(s) 
represents  the  feed  back  transfer  function. 

The  transient  response  of  the  system  is  characterized  by  the 
poles  of  Y(s)— which  are  the  roots  of  the  characteristic  equation 

1 + G(s)H(s)  = 0 (4-6) 

or,  in  our  case, 

1 + (-  ZPR  ttM  (G22  G24  + G21  G23))  = 0.  (4-7) 

Linear  stability  is  assured  if  all  the  poles  of  the  system  transfer 
function  are  in  the  left  half  plane,  i.e.,  if  roots  of  the  characteristic 
equation  are  negative  or  have  negative  real  parts. 

Using  the  program  CC,  the  roots  of  the  characteristic  equation, 
Eq.  4,7,  are  obtained  and  are  given  below: 

Zeros  of  1 - ZPR  Qm  (G22  G24  + G21  G23)  = 0 

- 2.1850  E-4 


- 7.8550  E-2 
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- 1.830 

+ 0.5463  ± j 4.0518 

- 50.4184 

-5.3712  ±j  47.4075 

- 3.6892  ± j 54.2463 

- 338.486 

- 26.888  ±J  13 1.7 18  and, 

- 900.0180 

The  above  values  for  the  roots  of  the  characteristic  equation  are 
obtained  for  the  following  BPGCR  operating  conditions: 

Core  Thermal  Power  = 100  MWt, 

Core  Average  Fuel  Temperature  = 1500  K, 

Average  Beryllium  Temperature  = 800  K. 

Average  Radiator  Temperature  = 1200  K, 

Core  Inlet  Stagnation  Pressure  = 31  atm. 

Fuel  Mass  Flow  Rate  = 140  Kg/s, 

Fuel  Mass  Feedback  Coefficient  = 0.035  Ak/k  per  Kg, 

Fuel  External  Loop  Circulation  Time  = 0.1  sec,  and 
Fuel  Core  Residence  Time  = 0.08  sec 
It  can  be  immediately  seen  that,  because  of  the  positive  root  (complex 
conjugate  pair)  on  the  right  half  plane,  the  system  is  unstable,  i.e,,  as 
time  progresses  after  a reactivity  perturbation,  the  disturbance 
(system  output)  grows  without  limit. 


162 


Earlier  in  this  chapter,  when  the  nonlinear  analysis  of  the 
system  was  performed  using  DYNAM,  it  was  found  that  the  fuel  mass 
feedback  is  a strong  feedback  mechanism,  and  the  system  can  be 
rendered  stable  or  unstable  depending  on  the  value  of  this  feedback 
coefficient.  It  is  decided  to  see  whether  the  fuel  mass  feedback  has 
the  same  effect  in  a linearized  model  also. 

To  verify  this,  the  value  of  the  fuel  mass  coefficient  is  reduced 
from  its  typical  value  of  0.035  to  a value  of  0.007  (a  factor  of  two 
reduction).  The  roots  of  the  characteristic  equation  are  again 
determined  using  CC  and  are  given  below; 

- 2.1850  E-4 

- 7.2280  E-2 

- 1.830 

- 0.4687  ±j  3.336 

- 52.177 

- 4.4060  ±j  26.2160 

-3.6810  ±j  54.2480 

- 338.486 

- 26.888  ± j 131.717  and, 

- 898.1870 

It  can  be  seen  that  the  system  is  now  stable  for  this  value  of  the  fuel 
mass  feedback  coefficient. 

It  is  sometimes  much  easier  to  visualize  the  system  stability  by 
looking  at  the  system  response  in  time  domain.  Using  CC,  the  time 
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response  of  the  BPGCR  (linearized  model)  is  studied  for  various  values 
of  the  fuel  mass  feedback  coefficient,  and  are  given  in  Figs.  4.30 
through  4.33.  Since  these  results  are  obtained  using  a linearized 
model  only  a small  reactivity  step  (10  cents)  is  applied  (large  reactivity 
perturbations  will  render  the  linearized  model  invalid).  Comparing 
these  figures  with  Figs.  4.10  through  4.12,  it  can  be  seen  that  the 
linearized  model  predicts  the  same  type  of  BPGCR  response  as  the 
nonlinear  model,  i.e.,  the  system  becomes  unstable  for  both  very  large 
and  very  small  values  of  the  fuel  mass  feedback  coefficient.  There  is  a 
range  of  values  of  this  coefficient  where  the  system  is  stable  to 
external  reactivity  perturbations.  The  actual  values  of  the  fuel  mass 
feedback  coefficient  (both  high  and  low)  at  which  the  system  becomes 
unstable  as  predicted  by  the  two  models  (linear  and  nonlinear), 
however,  are  different. 

Now  that  it  is  determined  that  the  system  is  stable— at  least  for 
certain  regimes  of  the  fuel  mass  feedback  coefficient,  the  next  step  is 
to  determine  whether  this  is  true  for  all  values  of  the  system  gain 
(reactor  power).  This  can  be  determined  from  a root  locus  plot  of  the 
transfer  function  of  the  system.  The  root  locus  method  is  based  on  the 
relationship  between  the  poles  and  zeros  of  the  closed  loop  transfer 
function  and  those  of  the  open  loop  transfer  function.  The  quantity 
’G(s)H(s)'  in  Equation  4.6  is  called  the  open  loop  transfer  function,  and 
it  plays  a very  important  part  in  the  stability  analysis  of  the  system  (see 
Ref.  26  for  details  of  the  root  locus  method).  Basically,  in  the  root 
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Figure  4-30.  Linear  Analysis  Prediction  of  the  BPGCR  Thermal  Power  Vs. 
Time  After  an  External  Reactivity  Insertion  - Case  I 
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Figure  4-31.  Linear  Analysis  Prediction  of  the  BPGCR  Thermal  Power  Vs. 
Time  After  an  External  Reactivity  Insertion  - Case  II 
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Figure  4-32.  Linear  Analysis  Prediction  of  the  BPGCR  Thermal  Power  Vs. 
Time  After  an  External  Reactivity  Insertion  - Case  III 
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Figure  4-33.  Linear  Analysis  Prediction  of  the  BPGCR  Thermal  Power  Vs. 
Time  after  an  External  Reactivity  Insertion  - Case  IV 
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locus  method,  the  behavior  of  the  roots  of  the  characteristic  equation 
(the  loci)  are  investigated  as  a function  of  the  system  gain.  A necessary 
and  sufficient  condition  for  linear  stability  is  that  the  roots  all  have 
negative  real  parts.  If  we  write  the  characteristic  equation  in  the  form 
1 + K G(s)  H(s)  = 0. 

then  the  root  locus  method  gives  the  loci  of  its  zeros  as  K varies  from 
zero  to  infinity  (the  gain  K include  the  reactor  power  and  the  fuel 
mass  feedback  coefficient  as  a part  of  it). 

The  open  loop  transfer  function  can  be  expressed  as  a ratio  of 
two  factored  polynomials 

rfci  Rfci  - K (s-ri)  (s-r2)  . . . (s-rj 

(s-pi)  (s-p2)  . . . (s-pm)  (4-8) 

where  'ri'  are  called  the  zeros  and  'pi'  are  called  the  poles.  Fig.  4.34 

shows  the  root  locus  plot  of  the  overall  BPGCR  transfer  function  for 
the  same  operating  conditions  given  earlier.  The  'X's  in  the  figure 
denotes  the  poles  and  the  'O's  denotes  the  zeros  of  the  open  loop 
transfer  function.  There  are  three  asymptotes,  approaching  infinity 
along  180,  +60,  and  -60  degrees  (Ref.  26  gives  the  rules  for 
constructing  a root  locus  plot).  In  the  figure,  for  the  sake  of  clarity, 
the  real  and  imaginary  axes  are  not  drawn  to  scale.  It  is  meant  for  the 
purpose  of  illustrating  the  behavior  of  the  roots  of  the  characteristic 
equation  only. 

From  the  figure  we  can  see  that  four  loci  cross  over  to  the  right 
hand  side  (positive  X coordinate  side)  as  the  system  gain  increases 


Figure  Not  to  Scale 
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Figure  4-34.  Root  Locus  Plot  of  the  BPGCR  with  Fuel  Mass  Feedback 
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from  zero  to  infinity.  The  system  gain  at  the  cross  over  points  are  also 
shown  in  Fig.  4.34.  It  can  be  seen  from  the  figure  that  the  limiting 
value  of  the  system  gain  is  about  2.5,  i.e.,  the  system  is  unstable  for  all 
values  of  gain,  K above  2.5.  Also,  from  the  root  locus  plot,  as  the  system 
gain  increases,  we  can  expect  the  oscillatory  behavior  of  the  system 
output  to  increase. 

So,  a linear  stability  analysis  also  indicates  that  the  BPGCR 
system  analyzed  in  this  dissertation  is  a stable  configuration  for  a 
range  of  operating  conditions  of  the  fuel  mass  feedback  coefficient. 


CHAPTER  5 

CONCLUSIONS,  RECOMMENDATIONS  AND 
AREAS  OF  FURTHER  RESEARCH 

Introduction 

The  d3niamic  analysis  of  a conceptual  gas  core  reactor  performed 
in  this  dissertation  indicates  that  the  Burst  Power  Gas  Core  Reactor 
system  or,  the  BPGCR  system,  is  capable  of  rapid  startup  and  very 
stable  operation  at  high  power  levels  in  a space  environment.  In  this 
dissertation,  the  BPGCR  is  primarily  considered  for  operation  in  a 
space  environment,  even  though  it  is  equally  capable  of  terrestrial 
operation. 

It  is  found  that  the  combination  of  the  gaseous  nature  of  the  fuel, 
plus  the  circulating  fuel  configuration  introduces  unique  feedback 
effects  which  tend  to  alter  significantly  the  dynamic  response  of  the 
BPGCR  compared  to  that  of  conventional,  solid  core  systems.  These 
unique  feedback  effects  are  also  found  to  make  the  system  very  stable 
during  operation  at  high  power  levels.  The  BPGCR  considered  in  this 
dissertation  is  designed  as  an  "open  flow"  system  with  nozzles,  rather 
than  valves  or  orifices,  at  the  core  inlet  and  at  the  core  exit.  In 
addition  to  the  simplicity  of  design,  such  a configuration  can  offer 
some  added  benefits  from  a materials  point  of  view  when  the  BPGCR  is 
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operated  using  highly  corrosive  gaseous  fuel  forms  like  UFe  or  UF4. 

The  analysis  performed  in  this  dissertation  indicates  that  this  BPGCR 
design,  with  only  nozzles  at  the  inlet  and  exit,  is  an  inherently  stable 
configuration  during  high  power  mode  of  operation. 

Results 

Transition  Analysis 

Since  the  BPGCR  is  a novel  system,  no  previous  analysis  exists  on 
the  dynamics  of  such  a system,  nor  is  any  information  available  on  the 
system  behavior  during  rapid  startups.  In  this  dissertation,  a computer 
program  (DYNAM)  is  developed  to  perform  this  analysis. 

The  computer  program  DYNAM  solves  simultaneously  the 
coupled  neutronics,  heat  transfer,  and  thermodynamic  equations  of 
the  BPGCR  attached  to  a balance-of-plant.  Since  no  similar  program 
exists,  various  sections  of  DYNAM  are  benchmarked  against  standard 
computer  codes,  or,  in  some  cases,  against  analytical  results.  DYNAM 
can  be  a very  valuable  tool  for  analyzing  the  startup  of  the  BPGCR  and 
also  for  further  refining  the  system  design. 

Using  DYNAM.  it  is  established  that  the  BPGCR  is  capable  of 
going  safely  and  in  a predictable  manner  from  a low  power  mode  of 
operation  to  a high  power  mode  (a  few  100s  of  MWt)  of  operation  in  a 
short  period  of  time  (of  the  order  of  100  seconds).  Additionally,  it  is 
found  that  only  simple  external  controlling  action  is  needed  to  control 
the  system  during  this  rapid  transition.  In  this  dissertation,  the 
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method  considered  is  to  control  the  system  by  controlling  the  core 
inlet  pressure. 

Different  startup  scenarios  (different  types  of  reactivity 
insertions)  are  considered  in  this  dissertation,  DYNAM  can  predict  the 
behavior  of  important  system  parameters  (like  core  pressure, 
temperature,  fuel  mass  flow  rates,  etc.)  during  any  of  these  startup 
scenarios.  An  advance  knowledge  of  the  behavior  of  these  parameters 
can  then  enable  one  to  select  desirable  transition  paths,  i.e.,  paths 
along  which  the  crucial  parameters  of  the  BPGCR  are  within 
satisfactory  bounds.  DYNAM  can  be  an  integral  part  of  a total  reactor 
control  system,  where  it  can  be  used  to  continuously  monitor  and 
adjust  the  inlet  core  pressure  so  as  to  implement  some  predeter- 
mined sequence  of  reactivity  insertion. 

The  computer  program  DYNAM  can  also  be  used  to  modify 
system  design  (for  example,  reactor  dimensions)  and  to  observe  the 
effect  of  such  modifications  on  important  system  variables  either 
during  transient  conditions,  or  during  normal  steady  state  operation. 
BPGCR  Stability  Analysis 

The  BPGCR  is  a circulating  fuel  reactor  using  gaseous  UFe  as  the 

fuel.  A circulating  fuel  reactor  by  itself  (i.e,,  with  no  feedback  effects)  is 
an  unstable  configuration.  This  means  that  any  perturbation  acting  on 
the  reactor  will  make  the  power  grow  without  limit  as  time 
progresses.  It  is  the  feedback  effects— which  are  always  present  in  a 
reactor  operating  at  power,  that  make  the  system  stable.  Because  of 
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the  gaseous  nature  of  the  fuel,  plus  the  circulating  fuel  configuration, 
the  BPGCR  possesses  some  unique  feedback  effects  which  are  found  to 
render  the  system  much  more  stable  than  conventional,  solid  core 
systems. 

To  accurately  predict  the  stability  of  the  BPGCR,  a nonlinear 
stability  analysis  of  the  system  is  performed  using  DYNAM,  Some  of  the 
more  important  results  obtained  from  this  analysis  are  given  below; 

a)  the  BPGCR  is  found  to  be  an  extremely  stable  configuration. 
Even  after  significant  reactivity  perturbations,  the  system  returns  to 
steady  state  within  a short  period  of  time.  It  is  observed  that  after  an 
external  reactivity  perturbation  of  the  order  of  $1.  the  resulting  power 
oscillations  have  amplitudes  of  no  more  than  30%  of  the  steady  state 
value. 

b)  the  inherent  BPGCR  response  to  external  reactivity 
perturbations  is  found  to  be  very  different  from  that  of  a conventional, 
solid  core,  non-circulating  fuel  reactor.  It  is  found  that  at  full  power, 
steady  state  conditions,  the  BPGCR  power  cannot  be  increased  by 
merely  applying  an  external  positive  reactivity  step  insertion.  To 
increase  the  BPGCR  power,  it  is  necessary  to  provide  additional 
positive  reactivity  during  the  transient  by  varying  a system  operating 
parameter  (e.g.,  by  increasing  the  core  inlet  pressure). 

c)  variations  in  the  external  loop  circulation  time  are  found  to 
influence  the  dynamic  behavior  of  the  BPGCR  for  the  examined 


175 


transients  and  the  combination  of  core  residence  time  and  external 
loop  circulation  times  considered  in  this  dissertation. 

d)  the  most  important  feedback  effect  that  influences  the 
stability  of  the  BPGCR  is  found  to  be  the  fuel  mass  (density)  effect. 
Depending  on  the  value  of  this  fuel  mass  feedback  coefficient,  the 
BPGCR  can  be  made  stable,  conditionally  stable  or  unstable. 

e)  limits  on  desired  values  of  the  fuel  mass  feedback  coefficient 
are  determined.  It  is  found  that  for  a certain  range  of  values  of  this 
coefficient,  the  reactor  exhibits  oscillatory  response  to  reactivity 
perturbations.  But  these  oscillations  are  heavily  damped  and,  by 
properly  selecting  the  operating  regime  of  the  BPGCR,  these 
oscillations  can  be  fully  eliminated— if  such  oscillations  are  deemed  to 
be  undesirable.  If  the  BPGCR  system  is  needed  to  provide  pulsed  mode 
power,  controlled  oscillations  in  the  power  output  can  help  to  reduce, 
or  even  eliminate,  some  of  the  power  conditioning  needs.  In  such  a 
situation,  the  value  of  the  fuel  mass  feedback  coefficient  can  be 
carefully  selected  to  provide  sustained  oscillatory  output. 

Sometimes  performing  a nonlinear  analysis  can  be  tedious,  and 
it  is  preferable  to  do  a linear  stability  analysis  if  such  an  analysis  can 
properly  predict  the  stability  of  the  system.  The  nonlinear  system 
equations  of  the  BPGCR  are  linearized  and  the  stability  of  the  BPGCR  is 
determined  using  the  established  methods  of  linear  stability  analysis. 

Comparing  the  results  of  the  linear  analysis  with  that  of  the 
nonlinear  analysis,  it  is  established  that  for  the  most  part,  a linear 
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stability  analysis  is  able  to  qualitatively  predict  correctly  the  stability  of 
the  BPGCR.  So,  for  initial  system  characterization,  a linear  stability 
analysis  may  prove  useful.  For  the  final  design  calculations,  however, 
one  needs  to  perform  a full-fledged,  nonlinear  analysis  to  obtain 
accurate  bounds  on  stability. 

Recommendations  for  Future  Work 

Foundations  for  the  dynamic  analysis  of  the  BPGCR  are  laid  in 
this  dissertation.  Some  refinements— which  will  make  the  analysis 
more  accurate,  are  discussed  below.  Some  of  the  recommendations 
are  short  term  and  others  need  a longer  time  to  implement.  On  a 
short  term  basis,  one  can  perform  the  foUoAving  tasks. 

a)  Some  further  parametric  studies  can  be  undertaken  to  study 
the  effect  of  changing  the  system  dimensions  or  design.  For  example, 
the  inlet  and  exit  throat  areas  can  be  varied,  and  the  effect  on  system 
dynamics  can  be  studied.  These  nozzle  areas  influence  the  fuel  mass 
inflow  and  outflow  rates  from  the  core  and,  hence,  will  influence  the 
reactivity  variations  during  transients.  The  nozzles  can  also  be 
replaced  by  valves  or  orifices  to  see  whether  they  will  make  the 
BPGCR  more  stable  or  not.  Similarly,  some  of  the  balance-of-plant 
variables  like  radiator  areas  can  be  varied  to  see  the  influence  of 
increasing  or  decreasing  the  external  loop  damping  on  the  d)mamic 
behavior  of  the  BPGCR. 

b)  In  order  to  reduce  computational  time,  only  one  delayed 
neutron  group  is  considered  in  DYNAM  while  solving  for  the 
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neutronics  equations.  A sample  case,  comparing  the  BPGCR  responses 
for  the  case  Avith  six  delayed  neutron  groups  versus  that  for  the  one 
delayed  neutron  group  are  given  in  Figs.  4.8  and  4.9.  These  figures 
indicate  that,  for  the  considered  transient,  the  results  from  the  two 
cases  are  almost  identical  and,  hence,  a one  delayed  neutron  group 
treatment  is  adequate.  For  certain  other  transients,  treatment  of  all  six 
delayed  neutron  groups  may  be  needed  to  accurately  predict  the 
reactor  power  behavior.  So,  for  final  design  calculations,  all  transients 
should  be  run  with  the  six  delayed  neutron  groups  option  of  DYNAM. 

c)  In  DYNAM,  the  external  loop  circulation  time  is  kept  fixed 
during  a transient.  It  is  not  allowed  to  change,  for  example,  as  a 
function  of  the  core  exit  velocity  of  the  fuel.  This  can  be  modified  by 
fixing  the  external  loop  length  and  letting  the  fuel  exit  velocity 
determine  the  external  loop  circulation  time.  This  might  influence  the 
d3mamic  response  of  the  BPGCR  by  influencing  the  core  inlet 
temperatures  and  delayed  neutron  precursors  reintroduced  into  the 
core. 

d)  Some  other  methods  of  controlling  the  BPGCR  can  also  be 
investigated.  In  this  dissertation,  the  BPGCR  is  controlled  by  means  of 
the  core  inlet  pressure.  How  the  system  will  respond  to  control  by 
means  of  the  fuel  mass  flow  rate  into  the  core  can  also  be  investigated. 
Control  of  the  BPGCR  by  means  of  core  exit  conditions  (e.g.,  pressure) 
can  also  be  studied,  although  this  might  be  more  difficult  to 


178 


implement  in  practice  because  of  the  extreme  conditions  (e.g.,  high 
fuel  temperature)  that  exist  at  the  core  exit. 

e)  In  the  balance -of-plant,  DYNAM  models  the  MHD  generator  as 
a turbine  generator.  This  can  be  refined  by  actually  treating  the  MHD 
generator  and  seeing  how  the  power  oscillations  inside  the  BPGCR 
affect  the  performance  of  the  MHD  generator  and  vice  versa. 

f)  It  is  found  in  this  dissertation  that  variations  in  the  fuel 
external  loop  circulation  time  influence  the  dynamic  behavior  (final 
steady  state  power,  amplitude  of  initial  power  oscillations,  etc.)  of  the 
BPGCR  after  an  external  reactivity  perturbation.  For  a circulating  fuel 
reactor,  it  is  actually  the  ratio  of  the  core  residence  time  and  the 
external  loop  circulation  time  that  is  important.  This  ratio  decides 
what  fraction  of  the  delayed  neutrons  is  lost  from  the  system  due  to 
the  circulation  of  delayed  neutron  precursors.  So,  other  values  for  the 
ratio  of  these  two  times  should  also  be  examined.  The  way  DYNAM  is 
written  currently,  there  is  only  provision  for  changing  the  external 
loop  circulation  time  of  the  fuel  (and  not  for  changing  the  core 
residence  time  of  the  fuel).  To  study  the  effect  of  varying  the  core 
residence  time  of  fuel  on  the  BPGCR  dynamics,  DYNAM  needs  to  be 
modified.  One  way  to  vary  the  core  residence  time  of  the  fuel  would  be 
to  artificially  vary  the  length  of  the  core  and,  hence,  the  transit  time  of 
the  fuel  through  the  core. 

Some  of  the  refinements  that  can  be  incorporated  into  the 
BPGCR  dynamic  analysis  on  a long  term  basis  are  discussed  next. 
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a)  In  DYNAM,  point  reactor  kinetics  (PRK)  equations  are  used  to 
predict  the  behavior  of  system  power  during  transients.  PRK  equations 
are  capable  of  predicting  the  behavior  of  the  system  power  with 
reasonable  accuracy,  but  a better  approach  would  be  to  use  space-time 
kinetics  to  predict  the  system  power  behavior.  In  DYNAM,  a uniform 
gas  density  is  assumed  throughout  the  core.  This  is  certainly  not  true. 
Unlike  in  a solid  core  system,  the  fuel  density  varies  significantly 
within  the  core  of  the  BPGCR. 

Since  the  fuel  is  flowing  continuously  through  the  core,  the  flow 
pattern  and  the  flow  rate  of  the  fuel  through  the  core  significantly 
influences  the  temperature  profile  and,  hence,  the  density  profile  of 
the  fuel  within  the  core.  Another  reason  for  the  non-uniform  density 
profile  within  the  core  is  the  fact  that  the  thermal  flux  peaks  near  the 
beryllium  moderator.  This  can  result  in  an  increased  power 
production  and,  hence,  increased  heat  generation  near  the  walls  of  the 
BPGCR  depending  upon  the  fuel  gas  density  distribution  within  the 
core.  Note  that  an  increased  heat  generation  near  the  walls  implies  a 
lower  gas  density  near  the  walls  which,  in  turn,  implies  a reduced 
fission  heat  generation  compared  to  the  center  of  the  core.  In 
practice,  this  effect  will  be  partly  cancelled  by  any  wall  cooling  that 
might  be  present,  i.e.,  some  coolant  will  be  flowing  through  the  walls 
to  keep  its  temperature  from  exceeding  the  melting  point  of  beryllium 
(or  BeO).  This  will  have  the  effect  of  making  the  fuel  denser  near  the 
wall.  It  is  apparent  from  the  above  discussion  that  the  gas  density 
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distribution  within  the  core  can  affect  the  spatial  flux  distribution 
within  the  core  and  vice  versa. 

Another  space-time  effect  that  might  be  present  in  the  system  is 
the  variation  of  the  core  average  flux  relative  to  that  of  the  reflector 
average  flux  during  a transient.  For  a cavity  configuration,  like  the 
BPGCR,  the  fission  neutrons  produced  within  the  core  get  moderated 
in  the  reflector  surrounding  the  core  (no  moderation  takes  place 
within  the  core  proper).  The  mean  time  taken  by  the  fast  neutrons  to 
return  to  the  core,  after  slowing  down  in  the  reflector,  is  only  of  the 
order  of  a few  milliseconds.  But,  when  considering  fast  transients,  the 
variation  of  the  core  average  flux  with  respect  to  that  of  the  reflector 
average  flux  needs  to  be  considered.  Such  space-time  effects  cannot  be 
properly  handled  by  PRK  models. 

b)  Another  phenomenon  that  needs  to  be  investigated  is  the 
effect  of  acoustic  oscillations  Avithin  the  BPGCR  core.  In  a cavity 
configuration  like  that  of  the  BPGCR,  employing  compressible  fluids 
and  with  pressure  dependent  internal  heat  generation,  it  is  an 
observed  fact  that  acoustic  oscillations  are  frequently  generated.  Such 
acoustic  oscillations  in  the  form  of  standing  pressure  waves  have  been 
predicted  for  gaseous  core  reactors.  These  acoustic  oscillations  can 
interact  with  the  neutronic  (amplitude)  oscillations  analyzed  in  this 
dissertaUon,  or  with  any  space-time  flux  oscillations  that  might  be 
present  in  the  system  (space-time  flux  oscillations  are  normally 


181 


observed  only  in  loosely-coupled  cores).  This  can  be  a significant 
effect,  and  needs  to  be  investigated. 

c)  A third  refinement  in  the  analysis  would  be  to  include  the 
second  half  of  the  BPGCR  (in  this  dissertation,  only  one  half  of  the 
BPGCR  is  considered),  and  introduce  the  coupling  effects  between  the 
two  halves  while  performing  the  dynamic  analysis.  Coupling  between 
two  such  cores  can  be  significant,  and  can  affect  the  d5niamic  response 
of  both  halves.  Panicker  (Ref.3)  has  shown  that,  for  the  pulsed  gas  core 
reactors  (PGCRs),  the  coupling  between  these  small  cores  can 
significantly  influence  their  dynamic  response. 


APPENDIX  A 

THE  COMPUTER  PROGRAM  DYNAM 
Introduction 

For  analyzing  the  transition  of  the  BPGCR  from  a low  power  mode 
of  operation  to  a high  power  mode  of  operation,  a computer  program 
called  DYNAM  is  developed  in  this  dissertation.  In  Chapter  3,  a 
detailed  description  of  this  program  and  the  models  used  in  it  are 
given.  DYNAM  is  also  used  to  perform  a nonlinear  stability  analysis  of 
the  BPGCR  in  Chapter  4.  This  Appendix  contains  a listing  of  the 
FORTRAN  program.  For  a detailed  description  of  the  models  used  in 
this  program,  refer  to  Chapter  3 and  Appendix  B of  this  dissertation. 
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PROGRAM  DYNAM 

INTEGER  N, KOUNT, K1 , KMAX, KOUNTl , KOUNT2 , NERR, KMAXl 
REAL*8  ATOMS (5) ,NUMDEN(4) ,DENSIT{5) , MASS (14) , H (7) , SECOND  (8) , 
1 KEFF , ENRICH, DEN ( 7 ) , POWER, DENEND ( 7 ) , MEANDN, CINLET ( 6 ) , 

1 COUT  (6,5500), MEANCP, MASSl , MASS2 , AVGMAS,  MASSIN,  MLEFT, 

1 TOTMAS,NHE,NUF6,  THROAT,  VOL,  HEIGHT,  DIA,  RADIUS,  AREAIN, 

1 AREAEX, AREAl , BETHIK, BEMASS, CPBE, RHOBE, RUNIV, RFUEL, 

1 SBCONS,CPBYCV,EMISS,LSTAR,BETA(6)  ,ALAMDA(6)  ,MACHIN, 

1 CORVEL, MACHSQ,  INRATE,  EXRATE,  RHOTOT,  RHCmX,  GAMMA, 

1 STAGPI, STAGPC, STAGTC, STAGTI, TINLET, TNEW, CORTMP, BETMP, 

1 COREPR, RATIO, FACTOR, FACTSQ, CYLINl , CYLIN2 , TERM  ( 3 ) , 

1 DELTAT, HEATNG (3) , STEP, TLOOP,  EXPO,  TIME,  TIME2,TIMAX,T2, 

1 T3,S0UNDI,DENIN,CRITPR,H2(7)  ,SECOND2(8) ,RHO(14) , 

1 RHOEXT, RHOINT, TEX, TRAD, T01,T04,TEXIT (5500) ,RSLOPE, 

1 EXITPR(5500) 

COMMON  /INTPRM/  LSTAR, BETA, ALAMDA 

CCayiMON  /CORDIA/  VOL,  HEIGHT,  DIA,  RADIUS,  AREAIN,  AREAEX,  THROAT, 

1 BETHIK, ENRICH 

COMMON  /CONST/  RUNIV, RFUEL, CPBYCV,  MEANCP,  CPBE,  SBCONS, EMISS 
COMMON  /TIMECO/  TIME, TIMAX, TLOOP, STEP 

C OPEN  (UNIT=1,FILE='FINAL2.0P1',STATUS='NEW) 

C OPEN  (UNIT=2,FILE=’FINAL2.0P2' ,STATUS=’NEW' ) 

C OPEN  (UNIT=3,FILE='FINAL2.0P3' ,STATUS='NEW' ) 

OPEN  (UNIT=4,FILE='FINAL2.0P4',STATUS='NEW ) 

DATA  ATOMS  /I . 2712, 0 . 2243, 8 . 9730, 19 . 8688, 0 . 1214/ 

DATA  MASS  /1 . 0,  2 . 0,  3 . 0,  4 , 0,  5 . 0,  6 . 0,  7 . 0,  8 , 0,  9 . 0,  9 . 2,  9 . 4, 

1 9.6,9.8,10.0/ 

DATA  KOUNT, KOUNT2,Kl, KMAX, N /O, 0, 0, 1000, 7/ 

DATA  DEN (1 ), CINLET  /O . 90E+8, 6*0 . 0/ 

DATA  STAGPI, STAGTI, CORTMP, TEX, TRAD, TOl, BETMP, CORVEL  /35.0, 

1 1200.0,1500.0,2000.0,1200.0, 900.0, 800.0,30.0/ 

DATA  RHOTOT, GAMMA, RHOMAX, TIME2/0 . 003, 0 . 02, 0 . 00710,0 . 0/ 

RHOBE  = 1800.0 
RSLOPE  =0.05 
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DO  10  1=1, 4,1 

NUMDEN(I)  = ATCa^Sd)  / (VOL  * l.OE+6) 
CONTINUE 

DENSIT(5)  = AT01S(5) 

DO  25  J=l,14,l 
DO  20  I = 1,4,1 

DENSIT(I)  = NUMDEN(I)  * MASS ( J) 
CONTINUE 

CALL  CRITIC (DENSIT,KEFF) 

RHO(J)  = 1.0  - (1.0  / KEFF) 
CONTINUE 

CALL  TABLE  (RHO,MASS,RHOTOT,MASSl) 


DO  30  I = 2,N,  1 

DEN (I)  = BETA(I-l) *DEN(1) / (LSTAR*ALAMDA(I-1) * (l.O-RHOTOT) ) 
CONTINUE 

CYLINl  = DIA  / 2.0 

CYLIN2  = CYLINl  + (BETHIK  / 100.0) 

BEMASS  = RHOBE  * 3.1416  * (CYLIN2  **  2 - CYLINl  **2)  * HEIGHT 

AREAl  = 3.1416  * DIA  * HEIGHT 

CONTINUE 

TIME  = TIME  + STEP 

IF  (TIME  .GE.  1.5)  THEN 
GO  TO  65 

ELSE  IF  (TIME  .GE.  TLOOP)  THEN 
RHOTOT  = 0.0035 
ENDIF 

CALL  RK4  (DEN, CINLET,N, RHOTOT, CORVEL,DENEND) 

MEANDN  = (DEN(l)  + DENEND(l))  / 2.0 
IF  (DENEND(l)  .GT.  l.OE+9)  THEN 

PRINT  *, 'POWER  EXCEEDED  1000  MWS' 

GO  TO  90 
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ENDIF 

NUF6  = MASSl  * 0.868  / 0.3495 
NHE  = MASSl  * 0.132  / 0.004 

COREPR  = (NUF6  + NHE)  * RUNIV  * CORTMP  * 9.872E-6  / VOL 

CALL  ACTUAL (STAGPI, STAGTI, COREPR, CPBYCV, TINLET,MACHIN, NERR) 

IF  (NERR  .GT.  0)  THEN 
KOUNT2  = KOUNT2  + 1 

PRINT*,  'AT',  TIME, 'SECS',  KOUNT2, ' CORE  INLET 
1 STAGNATION  PRESSURE (S)  LOWER  THAN  COREPR! ! ! ! ' 

ENDIF 

CORVEL  = 0.1  * MACHIN  * DSQRT (CPBYCV  * RFUEL  * TINLET) 

STAGPC  = STAGPI 

C CALL  FLORTE  (CORTMP, COREPR, EXRATE) 

EXRATE  = (0.6847  * 1.013E+5  * COREPR  * THROAT)  / (DSQRT (RFUEL 
1 * CORTMP) ) 

C CALL  REACT  (TIME2, ALAMDA, STEP, GAMMA, RHOMAX, RHOTOT) 

CALL  TABLE  (RHO, MASS, RHOTOT, MASS2) 

INRATE  = EXRATE  + ( (MASS2  - MASSl)  / STEP) 

CALL  CONTRL  (INRATE, COREPR, TINLET, AREAIN, STAGPI) 

CALL  TEMP  ( TINLET , INRATE , EXRATE , STEP , AREAl , BEMASS , DENEND ( 1 ) , 

1 MASS2, CORVEL, COREPR, TEX, CORTMP, BETMP, HEATNG) 

IF  (TIME  .GE.  TLOOP)  THEN 


DO  52  I = 1,6,1 

CINLET(I)  = COUT(I,l)  * DEXP (-ALAMDA (I ) * TLOOP) 
52  CONTINUE 


CALL  INTEMP  (TEXIT (1) , T04, TOl, TRAD, STAGTI, DELTAT) 
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IF  (STAGTI  .LE.  600.0)  THEN 
STAGTI  = 600.0 
ENDIF 

DO  54  I = 1,6,1 

DO  53  J = 1,  (KL-1) ,1 

GOUT  (I,  J)  = COUT(I,J+l) 

53  CONTINUE 

54  CONTINUE 

DO  55  I = 1,  (Kl-1)  ,1 
TEXIT(I)  = TEXIT(I+1) 

55  CONTINUE 

DO  56  I = 1,6,1 

COUT(I,Kl)  = DENEND(I+1) 

56  CONTINUE 
TEXIT(Kl)  = TEX 

ELSE 


K1  = K1  + 1 

DO  57  I = 1,6,1 

COUT(I,Kl)  = DENEND(I+1) 

57  CONTINUE 

TEXIT(Kl)  = TEX 

ENDIF 

KOUNT  = KOUNT  + 1 

IF  (KOUNT  .GE.  KMAX)  THEN 

C WRITE (1,92)  TIME,RHOTOT,MASS2,DENEND (1) 

C WRITE (2, 94)  TIME, COREPR, STAGPI,  CORTMP,  STAGTI, CORVEL, EXRATE 

C WRITE (3, 96)  TIME, TEX, T04, TOl, TRAD, STAGTI, DELTAT 

KOUNT  = 0 
ENDIF 

DO  60  I = 1,N 

DEN (I)  = DENEND(I) 
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60  CONTINUE 

MASSl  = MASS2 

GO  TO  50 

65  CONTINUE 

TIME  =0.0 

C DO  68  I = 1,K1,1 

C EXITPR(I)  = STAGPI 

C 68  CONTINUE 

C STAGPI  = STAGPI  +1.0 

RHOINT  = RHOTOT 

C WRITE(1,97)  TIME, RHOTOT, RHOINT,MASS2,DENEND(l) 

WRITE(4, 98)  TIME, TIME, RHOINT, DENEND (1) , STAGPI 
C WRITE (2,  94)  TIME, COREPR, STAGPI, CORTMP, STAGTI, CORVEL, EXRATE 

KMAX  = 200 

KOUNT  = 0 

KMAXl  = 500 

KOUNTl  = 0 

RHOEXT  = 0.00058 

RHOTOT  = RHOTOT  + RHOEXT 

70  CONTINUE 

TIME  = TIME  + STEP 
IF  (TIME  .GT.  TIMAX)  THEN 
GO  TO  90 
ENDIF 

CALL  RK4  (DEN, CINLET, N, RHOTOT, CORVEL, DENEND ) 

C MEANDN  = (DEN(l)  + DENEND(l))  / 2.0 

C IF  (DENEND (1)  .GT.  l.OE+9)  THEN 

C PRINT  *, 'POWER  EXCEEDED  1000  MWS' 
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C GO  TO  90 

C ENDIF 

NUF6  = MASSl  * 0.868  / 0.3495 
NHE  = MASSl  * 0.132  / 0.004 

COREPR  = (NUF6  + NHE)  * RUNIV  * CORTMP  * 9.872E-6  / VOL 
CALL  ACTUAL (STAGPI , STAGTI , COREPR, CPBYCV, TINLET, MACHIN, NERR) 

IF  (NERR  .GT.  0)  THEN 
KOUNT2  = KOUNT2  + 1 

PRINT*,  'AT',  TIME, 'SECS',  KOUNT2, ' CORE  INLET 
1 STAGNATION  PRESSURE (S)  LOWER  THAN  COREPR! ! ! ! ' 

ENDIF 

CORVEL  = 0.1  * MACHIN  * DSQRT (CPBYCV  * RFUEL  * TINLET) 
STAGPC  = STAGPI 

C CALL  FLORTE  (MACHSQ, FACTSQ, CORTMP, STAGPC, STAGTC, EXRATE) 

EXRATE  = (0.6847  * 1.013E+5  * COREPR  * THROAT)  / (DSQRT 
1 (RFUEL  * CORTMP) ) 

SOUNDI  = DSQRT (CPBYCV  * RFUEL  * TINLET) 

DENIN  = COREPR  * 1.013E+5  / (RFUEL  * TINLET) 

VELIN  = SOUNDI  * MACHIN 
CRITPR  = 0.528  * STAGPI 
IF  (COREPR  .LE.  CRITPR)  THEN 

INRATE  = (0.6847  * 1.013E+5  * STAGPI  * AREAIN)  / 

1 DSQRT (RFUEL  * STAGTI) 

ELSE 

INRATE  = DENIN  * AREAIN  * VELIN 
ENDIF 

MASS2  = MASSl  + STEP  * (INRATE  - EXRATE) 

CALL  TABLE  (MASS, RHO, MASS2, RHOINT) 

C RHOINT  = RHOINT  + RSLOPE  * (MASS2  - MASSl) 

CALL  TEMP  (TINLET, INRATE, EXRATE, STEP, AREAl, BEMASS, DENEND (1 ) , 
MASS2, CORVEL, COREPR, TEX, CORTMP, BETMP, HEATNG) 
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DO  72  I = 1,6,1 

CINLET(I)  = GOUT (1,1)  * DEXP ( -ALAMDA ( I ) * TLOOP) 

72  CONTINUE 

CALL  INTEMP  (TEXIT (1 ) , T04, TOl,  TRAD, STAGTI, DELTAT) 

DO  74  I = 1,6,1 

DO  73  J = 1, (Kl-l),l 

COUT(I,J)  = COUT(I,J+l) 

73  CONTINUE 

74  CONTINUE 

DO  75  I = 1, (Kl-1) ,1 
TEXIT (I)  = TEXIT (I+l) 

C EXITPR(I)=  EXITPR(I+1) 

75  CONTINUE 

DO  76  I = 1,6,1 

COUT(I,Kl)  = DENEND(I+1) 

76  CONTINUE 
TEXIT (Kl)  = TEX 

C EXITPR(Kl)  = COREPR 

C KOUNT  = KOUNT  + 1 

C IF  (KOUNT  .GE.  KMAX)  THEN 

C WRITE(1,97)  TIME,RHOTOT,RHOINT,MASS2,DENEND (1) 

C WRITE (2, 94)  TIME, COREPR, STAGP I,  CORTMP, STAGTI, CORVEL, EXRATE 

C WRITE (3, 96)  TIME,  TEX,  T04,  TOl,  TRAD, STAGTI, DELTAT 

C KOUNT  = 0 

C ENDIF 

KOUNTl  = KOUNTl  + 1 

IF  (KOUNTl  .GE.  KMAXl)  THEN 

WRITE (4, 98) TIME,  TIME, RHOINT,  DENEND (1) , STAGPI 
KOUNTl  = 0 
ENDIF 


DO  80  I = 1,N 
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DEN (I)  = DENEND(I) 

80  CONTINUE 

MASSl  = MASS2 
RHOTOT  = RHOINT  + RHOEXT 
C STAGPI  = EXITPR(l) 

GO  TO  70 

90  CONTINUE 

92  FORMAT(2F10.5,F10.3,2E12.4) 

94  FORMAT(F10.4,6F10.3) 

96  FORMAT(F10.4, 6F10.3) 

97  FORMAT(3F10.5,F10.3,E12.5) 

98  FORMAT(F10.3,1X,1H,  ,F10.5,1X,1H,  ,F10.5,1X,1H, ,E12.3,1X,1H, ,F8.2) 

END 

SUBROUTINE  CRITIC (DENSIT, KEFF) 

C THIS  SUBROUTINE  CALCULATES  THE  CORE  EFFECTIVE  MULTIPLICATION 
C FACTOR  FOR  A GIVEN  NUMBER  DENSITY  OF  FUEL  AND  BERYLLIUM  REFLECTOR. 

C IT  USES  A 2 GROUP,  2 REGION  DIFFUSION  THEORY  CALCULATION  ON  A 
C SPHERICAL  SYSTEM  FOR  ESTIMATING  THE  SYSTEM  KEFF.  THE  MICROSCOPIC 
C GROUP  CONSTANTS  - WHICH  THE  PROGRAM  USES  FOR  CALCULATING  THE  SYSTEM 
C KEFF  ARE  OBTAINED  VIA  INDEPENDENT  XSDRNPM  CALCULATIONS  AND  ARE  FED 
C TO  THE  PROGRAM  AS  INPUT.  SO, IF  THERE  ARE  ANY  DRASTIC  CHANGES  IN  THE 
C SYSTEM  GECMETRY  OR  FUEL  CONCENTRATION,  THESE  CONSTANTS  NEED  TO  BE 
C REEVALUATED . 

INTEGER  INTRVL(2) ,N1,N2,N 

REAL*8  SORCEl (400) ,SORCE2 (400) ,GP2SRC(400) ,GP1FLX(400) , 

1 GP2FLX (400) , DENSIT (5) ,TOTALl (5) ,TOTAL2 (5) ,ABS1(5) , 

1 ABS2 (5) ,FISS1 (5) ,FISS2 (5) , TRANSl (5)  ,TRANS2 (5) ,ONET01 (5) , 

1 ONET02 (5) ,TWOT02 (5) , SIGFl (2) , SIGF2 (2) , SIGRl (2) , SIGR2  (2) , 

1 SIGAl  (2) , SIGA2  (2) , SIGH  (2) , SIG12  (2) , SIG22  (2)  , SIGTRl  (2) , 

1 SIGTR2,D1 (2) ,D2 (2) , STEP (2) , NUl (2)  ,NU2 (2) , NU (2) , ENRICH, 

^ KOLD, KNEW, RADIUS, BETHIK,BTHIK2,  EPS,  TESTl,OLDSUM, 

1 NEWSUM,TERM(6) ,A1 (400) ,B1  (400) , Cl (400) , El (400) , FI (400) , 


1 

1 


A2 (400) ,B2 (400) ,C2 (400) , E2 (400)  ,F2  (400) , KEFF,  VOL, HEIGHT 
DIA, AREAIN, AREAEX, THROAT 
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1 

1 


CCMION  /MICRO/  T0TAL1,T0TAL2,FISS1,FISS2,ABS1,ABS2,NU1,NU2, 
TRANSl , TRANS2 , ONETOl , ONET02 , TWOT02 
COMMON  /MACRO/  Dl, D2, SIGAl, SIGA2, SIGFl, SIGF2, SIGH, SIG12, SIG22, 
SIGTRl, SIGTR2, SIGRl, SIGR2,  NU 

COMMON  /CORDIA/  VOL, HEIGHT, DIA, RADIUS, AREAIN, AREAEX, THROAT, 
BETHIK, ENRICH 


DATA  INTRVL,KOLD,EPS  /70, 200, 1 . 0, 1 . OE-4/ 


CALL  GRPCON  (DENSIT, ENRICH) 

BTHIK2  = BETHIK  + (0.71  / SIGTR2) 
STEP(l)  = RADIUS  / INTRVL(l) 

STEP (2)  = BTHIK2  / INTRVL(2) 

N1  = INTRVL(l) 

N2  = INTRVL(2) 

N = N1  + N2 


CALL  COEFF(N,Nl,  STEP,  RADIUS,  Dl  (1)  ,D2  (1)  , SIGRl  (1)  ,SIGR2  (1)  ,A1,B1, 
1 Cl, El, FI) 

CALL  COEFF  (N,N1,  STEP,  RADIUS,  Dl  (2)  ,D2  (2) , SIGRl  (2)  ,SIGR2(2)  ,A2,B2, 
1 C2,E2,F2) 

DO  110  I = 1,N1,1 
SORCEl(I)  =1.0 
110  CONTINUE 

DO  120  I = (Nl+1) , (N-1)  ,1 
SORCEl(I)  =0.0 
120  CONTINUE 

TERM(l)  = (4.0  /3.0)  * 3.1416  * ( STEP(l)  / 2.0)**3 

TERM(2)  = (3.1416/3.0)  * (STEP(1)**3) 

TERM(3)  = 3.1416  * ((2.0  * STEP(l)  * RADIUS  * RADIUS)  - 

(RADIUS  * STEP(l)  *STEP(D)  + (STEP(1)**3  / 6.0)) 
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OLDSUM  = SORCEl(l)  * TERM(l) 

DO  130  1=1, (Nl-l),l 

OLDSUM  = OLDSUM  + SORCEl(I)  * TERM(2)  * (1,0  + 12.0  * I * 
130  CONTINUE 

OLDSUM  = OLDSUM  + SORCEl(Nl)  * TERM (3) 

150  CONTINUE 

DO  155  I = 1, (N-1) ,1 

SORCEl(I)  = SORCEl(I)  / KOLD 
155  CONTINUE 

CALL  GRPFLX (N, N1 , SORCEl , A1 , B1 , Cl , El , FI , GPIFLX) 

DO  160  1=1, (Nl-l),l 

GP2SRC(I)  = SIG12(1)  * GPIFLX (I) 

160  CONTINUE 

GP2SRC(N1)  = (SIG12(1)  + SIG12(2))  / 2.0  * GPIFLX (Nl) 

DO  170  I = (Nl+1), (N-l),l 

GP2SRC(I)  = SIG12(2)  * GPIFLX (I) 

170  CONTINUE 

CALL  GRPFLX (N, Nl , GP2SRC, A2 , B2 , C2 , E2 , F2 , GP2FLX ) 

TERM(4)  = NU(1)  * SIGFl(l) 

TERM(5)  = NU(2)  * SIGFl (2) 

DO  175  I = 1,N1,1 

SORCE2(I)  = TERM(4)  * GPIFLX(I)  + TERM(5)  * GP2FLX(I) 
175  CONTINUE 

DO  180  I = (Nl+1) , (N-1) ,1 
SORCE2(I)  =0.0 
180  CONTINUE 

NEWSUM  = SORCE2(l)  * TERM(l) 

DO  185  I = 1, (Nl-1) ,1 

NEWSUM  = NEWSUM  + SORCE2(I)  * TERM(2)*(1.0  + 12.0  * I * I) 
185  CONTINUE 
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NEWSUM  = NEWSUM  + (S0RCE2(N1)  * TERM(3)  ) 

KNEW  = KOLD  * NEWSUM  / OLD SUM 
TESTl  = DABS ({KNEW  - KOLD)  / KNEW) 

IF  (TESTl  .LE.  EPS  ) THEN 
GO  TO  200 
ENDIF 

DO  190  I = 1,  (N-1)  ,1 
SORCEl(I)  = SORCE2(I) 

190  CONTINUE 

KOLD  = KNEW 
OLDSUM  = NEWSUM 
GO  TO  150 

200  CONTINUE 

KEFF  = KNEW 

RETURN 

END 

SUBROUTINE  GRPCON (DENSIT, ENRICH) 

C THIS  SUBROUTINE  CALCULATES  THE  MACROSCOPIC  GROUP  CONSTANTS  WHICH 
C ARE  NEEDED  BY  THE  CRITICALITY  PROGRAM  FOR  ESTIMATING  THE  SYSTEM 
C KEFF. 

REAL*8  ENRICH, DENSIT (5) ,TOTALl (5) ,TOTAL2 (5) ,ABS1 (5) ,ABS2(5) , 

1 TRANSl (5) ,TRANS2 (5) , FISSl (5) , FISS2 (5) ,ONET01 (5) , 

1 ONET02 (5) ,TWOT02 (5) , SIGTOT (2) , SIGRl (2) ,D1 (2) ,D2 (2) , 

1 SIGR2  (2) , SIGAl (2) , SIGA2 (2) , SIGFl (2) , SIGF2 (2) , SIGTRl (2) , 

1 SIGTR2,SIG11 (2)  ,SIG12(2)  ,SIG22 (2)  ,NU1 (2)  ,NU2 (2)  ,NU(2) 

COMMON  /MICRO/  TOTALl,  TOTAL2,  FISSl,  FISS2,  ABS1,ABS2, NUl, NU2, 

1 TRANSl , TRANS2 , ONETOl , ONET02 , TWOT02 

COMMON  /MACRO/  D1,D2, SIGAl, SIGA2, SIGFl, SIGF2, SIGH, SIG12, SIG22, 
1 SIGTRl, SIGTR2, SIGRl, SIGR2,NU 

DO  200  I = 1,2 
SIGAl (I)  =0.0 
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SIGFl(I)  =0.0 
SIGTR1(I)=  0.0 
200  CONTINUE 

SIGll(l)  = 0.0 
SIG12(1)  = 0.0 
SIG22(1)  =0.0 
SIGTOT(l)  =0.0 

DO  250  I = 1,4,1 

SIGTRl(l)  = SIGTRl(l)  + DENSIT(I)  * TRANSl(I) 
SIGTR1(2)  = SIGTRK2)  + DENSIT(I)  * TRANS2(I) 
SIGAl(l)  = SIGAl(l)  + DENSIT(I)  * ABSl(I) 

SIGA1(2)  = SIGA1(2)  + DENSIT(I)  * ABS2(I) 

SIGFl(l)  = SIGFl(l)  + DENSIT(I)  * FISSl(I) 

SIGF1(2)  = SIGF1(2)  + DENSIT(I)  * FISS2(I) 

SIGll(l)  = SIGll(l)  + DENSIT(I)  * ONETOl (I) 

SIG12(1)  = SIG12(1)  + DENSIT(I)  * ONET02(I) 

SIG22(1)  = SIG22(1)  + DENSIT(I)  * TWOT02(I) 

SIGTOT(l)  = SIGTOT(l)  + DENSIT(I)  * TOTALl(I) 
250  CONTINUE 

Dl(l)  = 1.0  / (3.0  * SIGTRl(D) 

Dl(2)  = 1.0  / (3.0  * SIGTR1(2)) 

D2(l)  = 1.0  / (3.0  * DENSIT(5)  * TRANSl (5) ) 

D2(2)  = 1.0  / (3.0  * DENSIT(5)  * TRANS2(5)) 

SIGA2(1)  = DENSIT(5)  * ABSl (5) 

SIGA2(2)  = DENSIT(5)  * ABS2(5) 

SIGF2(1)  =0.0 
SIGF2(2)  =0.0 

SIGTR2  = DENSIT(5)  * (TRANSl (5)  + TRANS2 (5) ) / 2.0 
SIGH  (2)  = DENSIT(5)  * ONETOl  (5) 

SIG12(2)  = DENSIT(5)  * ONET02 (5) 

SIG22(2)  = DENSIT(5)  * TW0T02(5) 

SIGTOT(2)  = DENSIT(5)  * TOTALl (5) 

SIGRl(l)  = SIGTOT(l)  - SIGll(l) 

SIGRK2)  = SIGA1(2) 

SIGR2(1)  = SIGTOT(2)  - SIGH  (2) 

SIGR2(2)  = SIGA2(2) 
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NU(1)  = ENRICH  * NUl(l)  + (1.0  - ENRICH)  * NU1(2) 

NU(2)  = ENRICH  * NU2(1)  + (1.0  - ENRICH)  * NU2(2) 

RETURN 

END 

SUBROUTINE  COEFF (N,N1, STEP, RADIUS, Dl, D2, SIGl, SIG2, A, B, C, E, F) 

C THIS  SUBROUTINE  CALCULATES  THE  ELEMENTS  OF  THE  COEFFICIENT  MATRIX 
C WHICH  IS  USED  FOR  SOLVING  THE  GROUP  FLUXES  (BOTH  FAST  AND  THERMAL) . 
C THIS  INFACT  RESULTS  IN  A PENTA  DIAGONAL  COEFFICIENT  MATRIX. 

INTEGER  N,N1 

REAL*8  STEP(2),A(400)  ,B(400)  ,C(400)  ,E(400) ,F(400) ,D1,D2,SIG1, 

1 SIG2,  DENOM,  FACT (7) , RADIUS, RI 

DO  300  I = 1,N 
A(I)  = 0.0 
B(I)  = 0.0 
C(I)  = 0.0 
E(I)  = 0.0 
F(I)  = 0.0 
300  CONTINUE 

FACT(l)  = Dl  / (2.0  * STEP(l)) 

FACT (2)  = D2  / (2.0  * STEP (2)) 

A(N1)  = -FACT(l) 

B(N1)  = 4.0  * FACT(l) 

C(N1)  =-3.0  * (FACT(l)  + FACT (2)) 

E(N1)  = 4.0  * FACT (2) 

F(N1)  = -FACT (2) 

FACT (3)  = Dl  / STEP(1)**2 
C(l)  =-FACT(3)  * 2.0  - SIGl 
E(l)  = FACT (3)  * 2.0 
DO  320  1=2, (Nl-1) ,1 
DENOM  =1.0/1 


B(I)  = FACTO)  * (1.0  - DENOM) 

E(I)  = FACTO)  * (1.0  + DENOM) 

C(I)  =-FACT(3)  * 2.0  - SIGl 
320  CONTINUE 

FACT (4)  = D2  / STEP (2) 

DO  340  I = (Nl+1) , (N-1) ,1 

RI  = RADIUS  + (I  - Nl)  * STEP  (2) 

FACT (5)  = (1.0  / STEP (2)  - 1.0  / RI) 

FACT (6)  = (1.0  / STEP (2)  + 1.0  / RI) 

B(I)  = FACT (4)  * FACT (5) 

E(I)  = FACT (4)  * FACT (6) 

C(I)  =-2.0  * FACT (4)  / STEP (2)  - SIG2 
340  CONTINUE 

DO  360  1=2, (N-2) ,1 
B(I)  = B(I)  / C(I-l) 

A(I+1)  = A(I+1)  / C(I-l) 

C(I)  = C(I)  - B(I)  * E(I-l) 

E(I)  = E(I)  - B(I)  * F(I-l) 

B(I+1)  = B(I+1)  - A(I+1)  * E(I-l) 

C(I+1)  = C(I+1)  - A(I+1)  * F(I-l) 

360  CONTINUE 

B(N-l)  = B(N-l)  / C(N-2) 

C(N-l)  = C(N-l)  - B(N-l)  * E(N-2) 

RETURN 

END 

SUBROUTINE  GRPFLX  (N, Nl,  R,  A,  B,  C,  E,  F,  X) 

C THIS  SUBROUTINE  SOLVES  FOR  THE  GROUP  FLUXES  USING  THE  COEFFICIENT 
C MATRIX  CALCULATED  BY  THE  PREVIOUS  SUBROUTINE.  THE  SOLUTION 
C TEC  NIQUE  IS  VERY  SIMILAR  TO  THE  THOMAS  ALGORITHM  - WHICH  IS  USED 
C FOR  SOLVING  TRIDIAGONAL  MATRICES.  IN  OUR  CASE  BECAUSE  OF  THE 
C SPECIAL  INTERFACE  CONDITON,  WE  END  UP  WITH  A PENTAD lAGONAL  MATRIX. 
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INTEGER  N,N1 

REAL*8  A (N-1 ) , B (N-1 ) , C (N-1 ) , E (N-1 ) , F (N-1 ) , R (N-1 ) , X (N-1 ) 

DO  400  1 = 1,  (N-1)  ,1 
R(I)  = -R(I) 

400  CONTINUE 

R(N1)  = 0.0 

DO  420  I = 2, (N-2) ,1 

R(I)  = R(I)  - B(I)  * R(I-l) 

R(I+1)  = R(I+1)  - A(I+1)  * R(I-l) 

420  CONTINUE 

R(N-l)  = R(N-l)  - B(N-l)  * R(N-2) 

X(N-l)  = R(N-l)  / C(N-l) 

X(N-2)  =(R(N-2)  - E(N-2)  * X(N-l))  / C(N-2) 

DO  450  I = (N-3),l,-l 

X(I)  = (R(I)  - X(I+1)  * E(I)  - X(I+2)  * F(D)  / C(I) 

450  CONTINUE 
RETURN 
END 

SUBROUTINE  SPLINE  (RHO, MASS, H, SECOND) 

C THIS  SUBROUTINE  SOLVES  FOR  THE  SECOND  DERIVATIVES  NEEDED  FOR  CALCU- 
C DATING  THE  COEFFICIENTS  OF  A CUBIC  SPLINE  FIT. 

REAL*8  RH0(8),MASS(8),H(7),TERM(7),A(6),B(6),C(6),R(6),G{6), 

1 SECOND (8) 

DO  500  I = 1,7,1 

H(I)  = RHO(I+l)  - RHO(I) 

TERM(I)  = (MASS(I+1)  - MASS (I))  / H(I) 

500  CONTINUE 


DO  520  I = 1,6,1 
A(I)  = H(I) 
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B(I)  = 2.0  * ( H(I+1)  + H(I)  ) 

C(I)  = H(I+1) 

R(I)  =(TERM(I+1)  - TERM(D)  * 6.0 
520  CONTINUE 

A(l)  = 0.0 

C(6)  = 0.0 

DO  540  I = 2,6,1 

A(I)  = A(I)  / B(I-l) 

B(I)  = B(I)  - A(I)  * C(I-l) 

R(I)  = R(I)  - A(I)  * R(I-l) 

540  CONTINUE 

G(6)  = R(6)  / B(6) 

DO  550  I = 5, 1,-1 

G(I)  = (R(I)  - G(I+1)  * C(I)  ) / B(I) 

550  CONTINUE 

SECOND (1)  =0.0 
SECOND (8)  =0.0 
DO  580  I = 1,6,1 

SECOND (I+l)  = G(I) 

580  CONTINUE 

RETURN 

END 

SUBROUTINE  TABLE (X, Y, IN, OUT) 

C THIS  SUBROUTINE  CALCULATES  THE  FUEL  MASS  IN  THE  CORE  CORRESPONDING 
C TO  A GIVEN  REACTIVITY  OR  THE  REACTIVITY  CORRESPONDING  TO  A GIVEN 
C CORE  FUEL  MASS  BY  LOOKING  UP  A TABLE.  A LINEAR  INTERPOLATION  IS  USED 
C FOR  INTERPOLATING  BETWEEN  THE  TABULATED  VALUES. 

INTEGER  I 

REAL*8  X(14) ,Y(14) ,IN,OUT 

IF  (IN  .LE.  X(2))  THEN 
1 = 1 

ELSE  IF  (IN  .LE.  X(3))  THEN 
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I 

= 2 

ELSE 

IF 

(IN  .LE.  X(4)) 

THEN 

I 

= 3 

ELSE 

IF 

(IN  .LE.  X(5)) 

THEN 

I 

= 4 

ELSE 

IF 

(IN  .LE.  X(6)) 

THEN 

I 

= 5 

ELSE 

IF 

(IN  .LE.  X(7)) 

THEN 

I 

= 6 

ELSE 

IF 

(IN  .LE.  X(8)) 

THEN 

I 

= 7 

ELSE 

IF 

(IN  .LE.  X(9)) 

THEN 

I 

= 8 

ELSE 

IF 

(IN  .LE.  X(10)) 

THEN 

I 

= 9 

ELSE 

IF 

(IN  .LE.  X(ll)) 

THEN 

I 

= 10 

ELSE 

IF 

(IN  .LE.  X(12)) 

THEN 

I 

= 11 

ELSE 

IF 

(IN  .LE.  X(13)) 

THEN 

I 

= 12 

ELSE 

I 

= 13 

ENDIF 

OUT  = Y(I)  + (((Y(I+1)  - Y(I)  )/(X(I+l)  -X(I)))  * (IN-X(I))) 

RETURN 

END 

SUBROUTINE  RK4 (BENIN, CIN, N, RHOTOT,VEL, DENOUT) 

C THIS  SUBROUTINE  PERFORMS  FOURTH  ORDER  RUNGE-KUTTA  INTEGRATION  WITH 
C CONSTANT  TIME  STEPS. 

INTEGER  N 

REAL*8  DENIN(7) ,CIN(6) ,DENOUT(7)  ,TEMP(7)  ,SLOPEl (7) ,SLOPE2 (7) , 
1 SLOPES (7) ,SLOPE4 (7) ,VEL, TIME, TIMAX, TLOOP, STEP, HFSTEP, 

1 LSTAR,BETA(6) ,ALAMDA(6) ,RHOTOT,  INVKEF,GENTME 

COMMON  /INTPRM/  LSTAR, BETA, ALAMDA 


COMMON  /TIMECO/  TIME, TIMAX, TLOOP, STEP 


INVKEF  = 1.0  - RHOTOT 

GENTME  = LSTAR  * INVKEF 

HFSTEP  = STEP  * 0.5 

CALL  RHS (BENIN, CIN,N, RHOTOT, GENTME, VEL, SLOPE!) 

DO  600  1=1, N 

TEMP(I)=DENIN(I)  + HFSTEP*SLOPEl (I) 

600  CONTINUE 

CALL  RHS (TEMP, CIN, N, RHOTOT, GENTME, VEL, SLOPE2 ) 

DO  620  1=1, N 

TEMP(I)=DENIN(I)  + HFSTEP*SLOPE2 (I ) 

620  CONTINUE 

CALL  RHS (TEMP, CIN,  N,  RHOTOT, GENTME, VEL, SLOPES) 

DO  640  1=1, N 

TEMP(I)=DENIN(I)  + STEP*SLOPE3 (I) 

640  CONTINUE 

CALL  RHS ( TEMP , CIN, N, RHOTOT , GENTME , VEL , SLOPE 4 ) 

DO  650  1=1, N 

DENOUT(I)=DENIN(I)+(STEP/6.)*(SLOPEl(I)  + 2.* 

1 SLOPE2(I)  + 2.*SLOPE3(I)  + SLOPE4(I)) 

650  CONTINUE 

RETURN 
END 

SUBROUTINE  RK42  (IN,T04) 

C THIS  SUBROUTINE  PERFORMS  FOURTH  ORDER  RUNGE-KUTTA  INTEGRATION  WITH 
C CONSTANT  TIME  STEPS. 

INTEGER  N 

REAL*8  IN(2) ,T04,TEMP(2) , SLOPE 1 (2) , SLOPE2 (2) , 

1 SLOPES (2) , SLOPE 4 (2) , TIME, TIMAX, TLOOP, STEP, HFSTEP 

COMMON  /TIMECO/  TIME, TIMAX, TLOOP, STEP 

N = 2 

HFSTEP  = STEP  * 0.5 

CALL  RHS2  (IN, N, TO 4, SLOPE!) 
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DO  660  I = 1,N 

TEMP  (I)  = IN(I)  + HFSTEP  * SLOPEl(I) 

660  CONTINUE 

CALL  RHS2  (TEMP,N,T04, SLOPE2) 

DO  670  I = 1,N 

TEMP  (I)  = IN(I)  + HFSTEP  * SLOPE2(I) 

670  CONTINUE 

CALL  RHS2  (TEMP, N, TO 4, SLOPES) 

DO  680  I = 1,N 

TEMP (I)  = IN(I)  + HFSTEP  * SLOPES (I) 

680  CONTINUE 

CALL  RHS2  (TEMP, N, T04, SLOPE4 ) 

DO  690  I = 1,N 

IN(I)  = IN(I)  + (STEP/6.0)  * (SLOPE! (I)  + 2.0  * 
1 SLOPE2(I)  + 2.0  * SLOPES (I)  + SLOPE4(I)) 

690  CONTINUE 

RETURN 
END 


SUBROUTINE  RHS  (DENIN, CIN, N, RHO, GENTME, VEL, DNDT) 

C THIS  SUBPROGRAM  COMPUTES  THE  R.H.S  OF  THE  SET  OF  NEUTRON  AND 
C PRECURSOR  DENSITY  EQUATIONS. 

INTEGER  N 

REAL*8  DENIN(N) ,CIN(6) ,DNDT(N) ,LSTAR,BETA(6) ,ALAMDA(6) , 

1 RHO, GENTME, VEL, TIME, TIMAX, TLOOP, STEP, FLOTRM(6) , SUM 

COMMON  /INTPRM/  LSTAR, BETA, ALAMDA 
COMMON  /TIMECO/  TIME, TIMAX, TLOOP, STEP 


SUM  =0.0 
DO  700  I = 1,6,1 

FLOTRM(I)  = (VEL  / 2.0)  * (CIN(I)  - DENIN(I+1)) 
700  CONTINUE 

DO  720  I = 1,6,1 

SUM  = SUM  + (ALAMDA (I)  * DENIN (I+l)) 


202 


720  CONTINUE 

DO  750  1=1, N,1 

IF  (I  .EQ.  1)  THEN 

DNDT(I)  = ( (RHO-0.0065)  * DENIN(I))  / GENTME  + SUM 
ELSE 

DNDT(I)  = (BETA(I-l)  * DENIN(l)  / GENTME)  - (DENIN(I)  * 

1 ALAMDA(I-D)  + FLOTRM(I-l) 

ENDIF 
750  CONTINUE 

RETURN 

END 

SUBROUTINE  RHS2  (IN, N, T04, DNDT) 

C THIS  SUBPROGRAM  COMPUTES  THE  R.H.S  OF  THE  SET  OF  HEAT  EXCHANGER  EXIT 
C TEMPERATURE  AND  AVERAGE  RADIATOR  TEMPERATURE  EQUATIONS. 

INTEGER  N,  I 

REAL*8  IN(2) ,DNDT(2) ,T04 

DO  800  1=1, N 

IF  (I  .EQ.  1)  THEN 

DNDT(I)  = -2.25  * IN(I)  + 2.50  * IN(I+1)  - 0.250  * T04 
ELSE 

DNDT(I)  = (0.3676  * IN(I-l))  - (5.88E-11  * (IN(I)**4)) 

1 - (0.7353  * IN(D)  + (0.3676  * T04) 

ENDIF 

800  CONTINUE 

RETURN 

END 

C SUBROUTINE  REACT (TIME2, ALAMDA, STEP, SLOPE, RHOMAX, RHO) 

C 

C THIS  SUBROUTINE  CALCULATES  THE  TOTAL  REACTIVITY  IN  THE  CORE  FOR 
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C EACH  TIME  STEP 
C 

C REAL*8  TIME2, EXPO, ALAMDA,RH0,RH01,RH02,RH0MAX, STEP, SLOPE 

C 

C IF  (TIME2  .GT.  0.0)  THEN 

C EXPO  = DEXP (- (ALAMDA  * TIME2) ) 

C RHOl  = RHOMAX  * EXPO 

C RH02  = 0.006458  * (1.0  - EXPO) 

C RHO  = RHOl  + RH02 

C TIME2  = TIME2  + STEP 

C ELSE 

C IF  (RHO  .GE.  0.0)  THEN 

C SLOPE  = 0.00015 

C ENDIF 

C RHO  = RHO  + (SLOPE  * STEP) 

C IF  (RHO  .GE.  RHOMAX)  THEN 

C TIME2  = TIME2  + STEP 

C ENDIF 

C ENDIF 

C 

C RETURN 

C END 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  FLORTE (MSQ, MACHEX, TEMP, STAGPC, STAGTC, MEX) 

THIS  SUBROUTINE  CALCULATES  THE  FUEL  MASS  OUTFLOW  RATE  FOR  A GIVEN 
CORE  STAGNATION  CONDITIONS  AND  BACK  PRESSURE. 


REAL*8  MSQ, MACHEX, MEX, PBACK, TEMP, STAGTC, STAGPC, 

1 TEXIT, VOL, AREAIN, AREAEX, THROAT, RUNIV, RFUEL, CPBYCV, 

1 BETHIK, ENRICH, HEIGHT, DIA, RADIUS, MEANCP, CPBE, RATIOl , 

1 RATIO, SBCONS,  EMISS 


COMMON  /CONST/  RUNIV, RFUEL, CPBYCV, MEANCP, CPBE, SBCONS, EMISS 
COMMON  /CORDIA/  VOL, HEIGHT, DIA, RADIUS, AREAIN, AREAEX, THROAT, 
BETHIK, ENRICH 


1 
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c 


c 


c 


MEX  = (0.6847  * 1.013E+5  * STAGPC  * THROAT) / (DSQRT (RFUEL 
* STAGTC) ) 


C 


1 


c 


RETURN 


C 


END 


SUBROUTINE  ACTUAL  (PO,  TO,  PI,  CPBYCV,  T1,MACH,NERR0R) 

C THIS  SUBROUTINE  CALCULATES  THE  ACTUAL  CONDITIONS  FOR  A GIVEN 
C STAGNATION  CONDITIONS . 

INTEGER  NERROR 

REAL*  8 XI , X2 , RATIO, MSQ, MACH, PO , TO , PI , T1 , CPBYCV 
XI  = CPBYCV  - 1.0 
X2  = XI  / CPBYCV 
RATIO  = PO  / PI 

MSQ  = ((RATIO**X2)  - 1.0)  * 2.0  / XI 

IF  (MSQ  .LT.  0.0)  THEN 
MSQ  = -1.0  * MSQ 
NERROR  = 1 

ELSE 

NERROR  = 0 
ENDIF 

MACH  = DSQRT (MSQ) 

T1  = TO  / (1.0  + (XI  * MSQ  / 2.0)) 

RETURN 

END 

SUBROUTINE  STAG (M, T1 , CPBYCV, TO ) 

C THIS  SUBROUTINE  CALCULATES  THE  STAGNATION  TEMPERATURE  GIVEN  THE 
C ACTUAL  CONDITIONS  AND  THE  SQUARE  OF  THE  MACH  NUMBER. 

REAL*8  M,T1, TO, CPBYCV, Y1,Y2, EXPO 
Y1  = CPBYCV  - 1.0 
Y2  = Y1  * M / 2.0 


EXPO  = CPBYCV  / Y1 
TO  = T1  * (1.0  + Y2) 

RETURN 

END 

SUBROUTINE  CONTRL (Ml , PRl , TEMP , AREA, PR2 ) 

C THIS  SUBROUTINE  CALCULATES  THE  INLET  STAGNATION  PRESSURE  NEEDED 
C FOR  OBTAINING  A GIVEN  MASS  INFLOW. 

REAL*8  M1,MACH,MSQ,  S1,S2,S3,S4,S5,PR1,PR2,RUNIV,RFUEL,MEANCP, 
1 CPBE, CPBYCV, SBCONS,EMISS, AREA, TEMP 

COMMON  /CONST/  RUNIV, RFUEL, CPBYCV, MEANCP , CPBE , SBCONS , EMISS 

51  = DSQRT (RFUEL  * TEMP  / CPBYCV) 

52  = Ml  / (AREA  * PRl  * 1.013E+5) 

MACH  = SI  * S2 

MSQ  = MACH* *2 

53  = CPBYCV  - 1.0 

54  = CPBYCV  / S3 

55  = 1.0  + (S3  * MSQ  / 2.0) 

PR2  = PRl  * (S5**S4) 

RETURN 

END 

SUBROUTINE  INTEMP  (TEX,  T04,  TOl,  TRAD, T02, DELTAT) 

C THIS  SUBROUTINE  CALCULATES  THE  CORE  INLET  STAGNATION  TEMPERATURE 
C GIVEN  THE  CORE  EXIT  TEMPERATURE  AND  BALANCE  OF  PLANT  CONSTANTS 
C LIKE  TURBINE  EFFICIENCY, COMPRESSOR  EFFICIENCY, TURBINE  ENTHALPY 
C EXTRACTION  RATIO  ETC. . 

INTEGER  N 

REAL*8  TEMPS (2) , TEX, TOl, T02, T04, ETAT, ETAC, EPSIL, UBARA,  TERM, 

1 DELTAT 

ETAT  =0.75 
EPSIL  =0.20 
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ETAC  =0.8 
UBARA  = 1.8E+5 

TERM  = 1 - (ETAT  / (ETAT  - EPSIL) ) 

TEMPS (1)  = TOl 
TEMPS (2)  = TRAD 

T04  = TEX  * (1  - EPSIL) 

CALL  RK42  (TEMPS, T04) 

TOl  = TEMPS (1) 

TRAD  = TEMPS (2) 

T02  = TOl  * (1  - TERM  / ETAC) 

DELTAT  = T02  - TOl 

RETURN 

END 

SUBROUTINE  TEMP (TIN, MIN, MEX, STEP, AREAl , BEMASS, POWER, MASS, 

1 VEL,COREPR,TEX,CORTMP,BETMP,HEATNG) 

C THIS  SUBROUTINE  CALCULATES  THE  CORE  FUEL  TEMPERATURE  AND  EXTERNAL 
C BERYLLIUM  TEMPERATURE  BASED  ON  POWER  GENERATED  WITHIN  THE  CORE 
C PLUS  RADIATIVE  AND  CONVECTIVE  HEAT  TRANSFER  FROM  THE  CORE  TO  THE 
C MODERATOR.  IT  ALSO  ACCOUNTS  FOR  THE  GAMMA  AND  NEUTRON  HEATING  OF 
C THE  MODERATOR. 

REAL*8  RENOLD, PRANDL, RUNIV, RFUEL,  CPBYCV, MEANCP, CPBE, SBCONS, 

1 EMISS, VOL, HEIGHT, DIA, RADIUS, AREAIN,AREAEX, THROAT, 

1 BETHIK, ENRICH, TIN, STEP, AREAl , BEMASS, POWER, MASS, MACHSQ, 

1 VEL , COREPR, CORTMP , BETMP , HEATNG ( 3 ) , TERM ( 7 ) , MOLWT , 

1 RHO,MU,CONDUT,MIN,MEX,TEX 

COMMON  /CONST/  RUNIV, RFUEL, CPBYCV, MEANCP, CPBE, SBCONS, EMISS 
COMMON  /CORDIA/  VOL, HEIGHT, DIA, RADIUS, AREAIN,AREAEX, THROAT, 

1 BETHIK, ENRICH 
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PRANDL  = 0.2355 
MOLWT  =28.2 

RHO  = COREPR  * MOLWT  / (0.08205  * CORTMP) 

MU  = 3.640E-8  * CORTMP  + 2.27E-5 
RENOLD  = RHO  * VEL  * DIA  / MU 
CONDUT  = 1.810E-4  * CORTMP  + 0.102 

TERM(l)  = POWER  * 0.9 

HEATNG(l)  = POWER  * 0.1 

TERM(2)  = AREAl  * SBCONS  * EMISS  * (CORTMP  *★  4 - BETMP  **  4) 
HEATNG(2)  = TERM (2) 

TERM(3)  = CONDUT  * 0.0364  * (RENOLD  * PRANDL  ) **  0.75  / DIA 
TERM(5)  = AREAl  * TERM(3)  * (CORTMP  - BETMP) 

HEATNG(3)  = TERM(5) 

TERM(6)  = HEATNG(l)  + HEATNG(2)  + HEATNG(3) 

TERM(7)  = MEANCP  * ( (MEX  * TEX)  - (MIN  * TIN)) 

TEX  = TEX  + (STEP  / (MASS  * MEANCP))  * (TERM(l)  - TERM(2) 

1 - TERM(5)  - TERM(7)) 

CORTMP  = (TEX  + TIN)  / 2 

C BETMP  = BETMP  + (STEP  / (BEMASS  * CPBE) ) * (TERM(6)  * 0.7) 

BETMP  = BETMP  + (STEP  / (BEMASS  * CPBE))  * (TERM(6)  * 1.0) 
RETURN 
END 

BLOCK  DATA 

REAL*8  TOTALl (5)  ,TOTAL2 (5) ,FISS1 (5)  ,FISS2 (5) ,ABS1 (5) ,ABS2  (5) , 
1 NUl (2) ,NU2 (2) ,TRANS1 (5)  ,TRANS2 (5)  ,ONET01 (5) ,ONET02 (5) , 

1 TWOT02 (5) ,D1 (2) ,D2 (2) ,SIGA1 (2) ,SIGA2 (2) , SIGFl (2) , 


1 


SIGF2  (2) , SIGH  (2) , SIG12  (2) , SIG22  (2) , SIGTRl  (2) , SIGTR2, 
SIGRl  (2),SIGR2(2)  ,NU(2),RHO(14)  ,MASS(14)  ,H(7)  ,SECOND(8) , 
LSTAR,BETA(6)  ,ALAMDA ( 6) , VOL,  AREAIN,AREAEX,  THROAT,  RADIUS, 
HEIGHT, DIA, BETHIK, ENRICH,  RUNIV,  RFUEL,  CPBYCV, MEANCP, 

TIME, TIMAX, TLOOP, STEP, SBCONS, CPBE, EMISS 


1 


1 


1 


1 


1 


COMMON  /MICRO/  TOTALl, TOTAL2, FISSl, FISS2,ABS1,ABS2, NUl, NU2, 
TRANSl , TRANS2 , ONETOl , ONET02 , TWOT02 
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COMMON  /MACRO/  Dl, D2,  SIGAl,  SIGA2,  SIGFl,  SIGF2, SIGH, SIG12, SIG22, 
1 SIGTR1,SIGTR2,SIGR1,SIGR2,NU 

COMMON  /INTPRM  / LSTAR, BETA, ALAMDA 

COMMON  /CORDIA/  VOL, HEIGHT, DIA, RADIUS, AREAIN, AREAEX, THROAT, 

1 BETHIK, ENRICH 

CCMMON  /CONST/  RUNIV, RFUEL, CPBYCV,MEANCP, CPBE, SBCONS, EMISS 
COMMON  /TIMECO/  TIME, TIMAX, TLOOP, STEP 

DATA  TOTALl,TOTAL2  /21 . 64, 23 . 39, 4 . 50, 1 . 804, 5 . 29, 

1 291.52,10.16,4.11,0.845,6.199/ 

DATA  FISS1,FISS2  /8 . 02, 8 . 27E-2, 3*0 . 0, 233 . 58, 4*0 . 0/ 

DATA  ABS1,ABS2  /12 . 31, 7 . 37, 5 . 02E-3, 5 . 682E-5, 4 . 125E-3, 

1 276.8,1.244,4.265E-3,3.107E-3,4.894E-3/ 

DATA  TRANS1,TRANS2  /7 . 901, 8 . 175, 3 . 173, 2 . 633, 3 . 716, 

1 291.52,10.16,3.971,0.2844,6.015/ 

DATA  ONETOl  /9 . 324, 16 . 02, 4 . 48, 1 . 794, 5 . 219/ 

DATA  ONET02  /3 . 36E-3, 2 . 376E-3, 1 . 331E-2, 1 , OlE-2, 7 . 89E-2/ 

DATA  TWOT02  /14 . 717, 8 . 917, 4 . 108, 0 . 842, 6 . 194/ 

DATA  NU1,NU2  /2 . 431, 2 . 781, 2 . 419, 2 . 320/ 

DATA  LSTAR  /6.0E-4/ 

DATA  BETA  /O . 000215, 0 . 00142, 0 . 00127, 0 . 00257, 0 . 00075, 0 . 00027/ 
DATA  ALAMDA  /O . 0124, 0 . 0305, 0 . Ill, 0 . 301 , 1 . 1, 3 . 0/ 

DATA  VOL,  HEIGHT,  DIA,  RADIUS,  AREAIN,  AREAEX,  THROAT,  BETHIK, 

1 ENRICH/1.57,2.0,1.0,72.1,0.09,0.09,0.05,50.0,0.85/ 

DATA  RUNIV, RFUEL, CPBYCV,MEANCP, CPBE, SBCONS, EMISS  /8.314, 

1 294. 98, 1.403, 1000. 0,2428. 46, 5. 67E-8, 0.25/ 

DATA  TIME, TIMAX, TLOOP, STEP  /O . 0, 10 . 0, 0 . 1, 0 . 0001/ 


END 


APPENDIX  B 

BPGCR  SYSTEM  EQUATIONS.  AND  THE  DERIVATION 
OF  THE  BPGCR  TRANSFER  FUNCTIONS 

Chapter  3 of  this  dissertation  describes  the  transition  analysis  of 
the  BPGCR  from  a low  power  mode  of  operation  to  a high  power  mode 
of  operation.  In  Chapter  3,  brief  descriptions  of  the  equations  used  for 
modelling  the  system  are  given.  This  Appendix  describes  the  full  set  of 
equations  as  used  in  Chapter  3 and  Chapter  4 of  this  dissertation. 
Chapter  4 includes  a linear  stability  analysis  of  the  BPGCR.  For  this 
analysis,  the  nonlinear  system  of  equations  are  linearized  and  transfer 
functions  are  derived  from  this  set  of  linearized  equations.  This 
Appendix  describes  the  procedures  used  for  linearizing  the  system 
equations,  and  gives  the  methodology  for  deriving  the  BPGCR  transfer 
functions.  It  also  contains  the  transfer  functions  used  in  Chapter  4.  A 
list  of  variables  used  in  this  Appendix  is  given  at  the  end  of  this 
Appendix. 
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Relationship  Between  Core  Inlet  Temperature.  Average 
Core  Temperature.  Core  Power  and  Core  Exit  Temperature 
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The  equations  relating  the  average  core  fuel  temperature  and 
the  average  beryllium  moderator  temperature  are  given  by 


Mf  (t)  = N(t)  - h Aw  (Tf  - Tbc)  - C^  (m3'(t)  Ts'  - m2  (t)  T2) 


(B-1) 


and. 


Mbc  C?"  = h Aw  (Tf  - Tee)  - Ci 

where  " Ci"  represents  a constant  heat  removal  term  from  the 


(B-2) 


beryllium  moderator,  i.e.,  it  is  assumed  that  a constant  amount  of  heat 
is  removed  from  the  beryllium  moderator/ reflector  at  all  times, 
irrespective  of  the  core  power  or  other  operating  conditions.  This 
could  be  accomplished,  for  example,  by  passing  a gas  such  as  helium 
through  the  beryllium  wall  at  variable  flow  rates.  All  other  variables 
mentioned  in  the  above  equations  are  described  at  the  end  of  this 
Appendix.  The  state  points  2'  and  3'  used  in  the  Equation  (B-1) 
represents  the  conditions  inside  the  BPGCR  core  just  after  the  inlet 
nozzle  and  just  before  the  exit  nozzle,  respectively.  Figure  B-1  shows 
the  state  points  (1,  2,  3,  etc.)  used  elsewhere  in  this  Appendix,  but  it 
does  not  show  the  points  2'  and  3'. 
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CO 


Figure  B- 1 . T-S  Diagram  for  t±ie  BPGCR  Power  System 
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And  finally,  the  core  inlet  temperature,  the  core  exit 
temperature  and  the  core  (fuel)  average  temperature  are  related  by 
the  expression 


Tf  = 


T2'  + T3' 
2 


(B-3) 


For  performing  a linear  stability  analysis  of  the  BPGCR,  we  have 
to  derive  the  transfer  functions  relating  the  variables  involved.  To 
obtain  the  transfer  functions,  the  above  set  of  equations  have  to  be 
linearized.  The  procedure  used  is  as  follows:  first  all  the  variables  in 
the  above  equations  are  written  as  a sum  of  two  parts— one  part  being 
their  steady  state  value  and  the  other  being  their  variations  from  their 
steady  state  values.  As  an  example,  the  core  thermal  power,  N(t),  is 
rewritten  as 

N(t)  = Nss  + 3N(t) 

where  'Nss'  denotes  the  steady  state  part  and  '9N'  denotes  the  time 
dependent  variation  from  the  steady  state  value  (in  a linear  stability 
analysis,  it  will  assumed  that  the  magnitude  of  the  time  dependent 
part  is  small  compared  to  the  steady  state  value).  All  other  variables 
are  similarly  written  as  a sum  of  a steady  state  value  and  a time 
dependent  variation  from  the  steady  state  value.  Equations  (B-1),  (B- 
2),  and  (B-3)  are  rewritten  in  terms  of  these  new  variables. 

Subtracting  the  steady  state  equations  from  this  new  set  of  equations 
results  in  a third  set  of  equations  given  below— which  involves  only  the 
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derivatives  of  the  time  dependent  part  of  all  the  variables  given  in 
Equations  (B-1).  (B-2),  and  (B-3). 


Mp  = 3N(t)  - h Aw  OTf  - 3TBe)  - 

[m«(3T3'  - BT2)  + Ti  am3'  - 1*2^  am2]  (B -4) 

Mse  C?"  = h Aw  OTf  - 9TBe)  (5.5) 


aTf  = 


3T2  + 9T3' 
2 


(B-6) 


In  arriving  at  Eq.  (B-4),  the  following  assumptions  have  been 

made; 


a)  The  fuel  mass  inside  the  core  is  taken  to  be  a constant,  M(ss, 
for  estimating  the  fuel  heat  capacity  term  (i.e.,  MfCpf  term  on  the  right 

hand  side  of  Eq,  (B-4)).  After  the  system  becomes  critical,  the  fuel 
mass  inside  the  core  varies  very  little  for  any  significant  reactivity 
variations  (see,  e.g..  Fig,  3.19).  This  is  due  to  the  fact  that  the  core 
reactivity  is  a strong  function  of  fuel  mass  inside  the  core  and,  hence, 
the  total  fuel  mass  inside  the  core  varies  only  slightly  from  its  steady 
state  value  for  the  kind  of  reactivity  variations  allowed  by  a linear 
stability  analysis.  (Note,  however,  that  it  has  not  been  assumed  that 
core  inlet  and  outlet  fuel  mass  flow  rates  are  the  same  or  they  remain 
constant  to  provide  a constant  fuel  mass  inside  the  core.) 

b)  Terms  involving  the  product  of  two  small  variations  from  the 
steady  state  values  (i.e.,  the  "9"  quantiies)  are  neglected  in  comparison 
to  other  terms. 
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Now  to  derive  the  transfer  functions  from  the  above  set  of 
equations,  the  Laplace  transform  of  Eqs.  (B-4),  (B-5)  and  (B-6)  are 
taken  to  obtain  Eqs.  (B-7),  (B-8)  and  (B-9)  given  below; 


_ P s 9Tf(s)  = =^i-^N(s)  - 9Tf  (s)  + aiBe(s)  - 
h h A\v 

^[m-OT3  (s)  - 3T2(s))  + Tf  8013  (5)  - 8012(8)] 

hAw 

(B-7) 

— s8TBe(s)  = 8ff(s)  - 8TBe(s) 

(B-8) 

hAw 

afrfs)  = 

2 * 

(B-9) 

These  equations  can  be  rewritten  as 

ai  s 3Tf{s)  = a2  3N(s)  - 9Tf  (s)  + 9TBe(s)  - as  - 3T2(s)]  - 


as  Ts^  8ois’  + as  l2^  8012 

(B-10) 

a4  s 8Tbc(s)  = 8Tf(s)  - 8Tbc(s) 

(B-11) 

afKs)  = 

2 

(B-12) 

where 

mP  ^ 

= -= ^ 

h A^  , 

a2  = — 

h A^  , 

as  = -- 

h A-^  , and 

and 


Equation  (B-11)  can  be  solved  for  8TBe  to  obtain 
^Be  ^ 1 

3Tf  a4S+l  . (B-13) 

(The  right  hand  side  of  Eq.  (B-13)  is  denoted  as  the  transfer  function 
G8  in  Chapter  4.  Similarly.  Eq.  (B-12)  provides  an  expression  for  the 
transfer  function  G7  used  in  Chapter  4,  viz.,  G7  = 1/2.)  Using  Eq.  (B- 
13)  to  eliminate  5TBe  from  Equation  (B-10)  and  using  the 

relationshiop  (B-12),  we  can  finally  obtain  the  relationship  between 
the  core  inlet  temperature,  the  core  exit  temperature,  the  core  inlet 
and  exit  fuel  mass  flow  rates  and  the  core  power  as 

as  = a68T2  + a20N  + a3T|®3m2  - asTf^dma  (B- 14) 

where 


^ ai  a4S  _ 

^5  = T ^ ^ 2(1  ^ a4  s)  ^ ^3”^ 


, and 


a«  = -^s- ^ 

^ 2 2(1  + a4  s) 


asmf 


Equation  (B-14)  can  be  rewritten  in  a form  compatible  with  the 
transfer  function  notation  used  in  Chapter  4 as 

^3'  ~ Gi0T2'  + G23N  + G39m2'  + G43m2'  (B-15) 


where 
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G2  = ^ 

^ as 


G3  = 


^ /X'SS 

as  T2 
as 


, and 


G4  = - 


as  T2' 
as 


Balance-of  Plant  Model  Equations 

In  the  BPGCR  system  being  analyzed  in  this  dissertation,  there  is 
a converging  nozzle  at  the  core  inlet  and  a converging-diverging  nozzle 
at  the  core  exit  (see  Fig.  3-3).  It  is  assumed  in  this  analysis  that  the 
flow  through  these  nozzles  is  isentropic  and,  hence,  there  is  no 
stagnation  temperature  (or  pressure)  change  in  going  through  these 
nozzles.  The  temperatures  T2'  and  Ta’  used  in  Equations  (B-1)  through 
(B-15)  represents  the  conditions  inside  the  core  Just  after  the  inlet 
nozzle  and  just  before  the  exit  nozzle,  respectively.  The  points  2 and  3 
shown  in  Fig.  B-1  represents  conditions  outside  the  core,  just  before 
the  inlet  nozzle  and  just  after  the  exit  nozzle,  respectively.  Because  of 
the  isentropic  flow  assumption,  we  can  write 


T2’  = T02'  = To2  , and 

3 = i03  = i03 


where  the  subscript  "0"  refers  to  stagnation  conditions.  When  the  fuel 
velocity  is  small  (a  few  tenths  of  a Mach  number),  the  total 
temperature  and  the  stagnation  temperature  can  be  interchanged  with 
very  little  error. 
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Turbine  Equations 

In  the  actual  BPGCR  system,  an  MHD  generator  is  employed  to 
convert  the  thermal  power  generated  Avithin  the  BPGCR  to  electrical 
power.  For  the  purpose  of  this  dissertation,  this  MHD  generator  is 
modelled  as  a conventional  turbine  (see  Chapter  3 for  more  discussion 
on  this).  So,  the  nozzle  and  diffuser— normally  associated  with  an  MHD 
generator,  are  not  included  here. 

For  a turbine,  the  pressure  ratio,  rt  is  given  by 

rt  = Po3/Po4  (B-16) 

(see  Fig.  B-1  for  the  location  of  the  state  points  3 and  4).  The  turbine 
enthalpy  extraction  ratio,  e,  and  isentropic  efficiency,  rit,  are  defined  as 


e = (To3  - To4)  / To3  . and  (B-1 7) 

■•It  = (Tq3  - T04)  / (Tq3  - To4i)  . (B-1 8) 


Using  Equation  (B-1 7)  and  the  isentropic  relationship  between 
pressure  and  temperature,  viz.. 


we  get 


Tq3  _ 
Tq41 


(B-19) 


To4  = Tq3  (1-  e),  and 
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JIlll 

To4t  = [rt]\  Y I To3  . 

Using  the  above  two  relations,  we  can  obtain  the  expression  for  the 
isentropic  turbine  efficiency  as 

^t  = 


[1-  m Y 


and,  hence. 


rt  = 


Tit 


Tit-e 


niJL 

Y-l 


So,  given  the  turbine  inlet  conditions  and  the  turbine  parameters  like 
isentropic  efficiency,  enthalpy  extraction  ratio,  etc.,  the  conditions  at 
the  turbine  exit  can  be  obtained  from  the  relations 


Tq4  = Tq3  (1  - e) 


(B-20) 


and 


Po4  = Po3 


Tit  -e 


(B-21) 

The  equation  for  the  turbine  exit  temperature,  Tq4,  (i.e.,  Eq.  B- 
20)  can  be  linearized  as  follows 
3To4  = (1-  e)  3Tq3  • 


Or,  using  the  Laplace  transform  variable. 


dTo4  (s)  = (1-  e)  dTTos  (s)  . 


(B-22) 
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Heat  Exchanger/Radiator  Equations 

After  exiting  from  the  turbine,  the  fuel/working  fluid  (gas 
mixture)  passes  through  a heat  exchanger  to  transfer  the  waste  heat  to 
a radiator.  It  is  assumed,  for  this  analysis,  that  an  average  heat 
exchanger  temperature  and  an  average  radiator  temperature  can  be 
used  to  describe  the  heat  transfer  between  the  two.  With  this 
assumption,  the  relevent  equations  for  the  average  heat  exchanger  and 
radiator  temperatures  are  as  given  below; 

[Mf^  cfp  + ^ cfp(T04  - To ,)  - 

t^hx  ^hx  ("I^hx  ■ ^rad) 

(B-23) 


[Msea,  ^etal  ^ = 


^hx  -^hx  ■ Trad)  " ^rad  £ Trad 


and. 


(B-24) 


Thx 


= Tqi  + Tq4 
2 


(B-25) 


As  was  done  for  the  core  fuel  temperature  equations  (Equations 
(B-1),  (B-2)  and  (B-3)  ),  the  Equations  (B-23),  (B-24)  and  (B-25)  can 
be  first  linearized  by  considering  small  variations  around  some  mean 
(i.e.,  rewriting  all  the  variables  as  a sum  of  a steady  state  term  and  a 
time  dependent  variation  around  this  steady  state  value)  and  then 
subtracting  the  steady  state  equation  from  this  new  equations.  Finally, 
the  Laplace  transform  of  the  equations  are  taken  to  obtain  the 
following  set  of  equations 
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[S]  5Th,  = (5to4  - 8To i)  - (sThx  - 5T,ad) 

STrad  (S)  (C4  +S)=  C3  5fhx  (S) 

5Thx  (s)  = + 5Tq4(s) 

2 

where 


(B-26) 

(B-27) 

(B-28) 


C3  = 


^hx  ^hx 


M^IFai  + Mhe  . and 

C4  = C2  + C5 


p ^»^ne  '-'p 

aT^ad' 


.aTradJss  with 


C5  = ^ ^ 

M^fFai  + Mhe  . 

(Note  that  before  Laplace  transform  of  Eq.  (B-24)  can  be  taken,  it  has 
to  be  linearized  because  of  the  nonlinear  term  involving  Trad"^-) 


The  three  Equations,  (B-26),  (B-27),  and  (B-28)  can  be  solved 
simultaneously  to  finally  obtain  an  expression  for  the  temperature  at 
the  outllet  of  the  heat  exchanger  (Toi)  in  terms  of  the  heat  exchanger 

inlet  temperature  (T04)  as  follows 


5Toi  = 


2 Ce  - C7  + 


2 C0  +C7  - 


C2  C7 
(C4  + s) 
C2C7 
(C4  + s) 


5To4 


(B-29) 


where 
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C2  . and 


Q_  _ hhx  Ahx 
C2 


or,  taking  the  Laplace  transform  of  Eq.  (B-29),  we  can  write  it  as 


Compressor  Equations 

The  isentropic  efficiency  of  the  compressor,  ric.  is  defined  as 

"He  = (Toi  - T021)  / (Toi  - T02} 

(see  Fig.  B-1  for  the  location  of  the  state  points  1 and  2) 
so. 


The  isentropic  relationship  between  pressure  and  temperature  is 
given  by 


The  pressure  ratio  on  the  right  hand  side  of  Eq.  (B-31)  is  known,  once 
the  turbine  parameters  (efficiency,  enthalpy  extraction  ratio,  etc.)  are 
given.  From  the  last  two  relations  given  above,  we  get 


6Toi(s)  = Cs  5To4(s)  . 


(B-30) 


r lY-l 

To2i  _ Pq2  ~T~ 


Toi  I.P01. 


(B-31) 


^02  ^ 1 . 


Y-1 


TOI  ^c 


Y 
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or,  in  terms  of  Laplace  transform  variable. 


dTp2(s) 

dToi(s) 


1-^ 

1- 

T-l 

Y 

^c 

L iPoii 

(B-32) 


Equations  (B-22),  (B-30)  and  (B-32)  can  be  combined  to  obtain 


the  complete  transfer  function  of  the  balance -of-plant,  except  for  the 


time  delay  functions,  G5,  in  Chapter  4 as 


STq2  (s) 
5To3(s) 


= (Cs)  (1  - e) 


Tit-e) 


(B-33) 


Time  Delay 

As  mentioned  in  Chapter  4,  the  total  delay  of  the  external  loop  is 
divided  into  two  parts;  one  for  the  transport  of  fuel  from  the  core  to 
the  heat  exchanger  and  the  other  for  the  transport  of  fuel  gas  mixture 
from  the  heat  exchanger  back  to  the  core. 

We  can  write  an  equation  for  the  time  delay  of  the  form 


To4(t)  = To3(t  - ti  ) 


where  'xi'  time  delay  for  the  transport  of  the  fuel  from  the  core  to  the 

heat  exchanger  and  T03'  and  T04'  are  the  core  exit  temperature  and 

the  heat  exchanger  inlet  temperature,  respectively. 

Taking  the  Laplace  transform  of  the  above  equation  gives 
5To4  (s)  = 3To3  exp(-xis) 

which  can  be  approximated  by  (for  small  values  of  xi)  as 
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3To4  (s)  = aTo3(s)  / (1  +Ti  s)  . (B-34) 

(The  typical  values  of  the  external  loop  circulation  times  used  in  this 
dissertation  are  of  the  order  of  a tenth  of  a second.  Since,  "xi " is  only 

half  of  this,  the  above  approximation  should  not  be  too  bad.) 

Similarly, 

3To2  (s)  = 3Toi  (s)  / (1  + X2  s) 

where  "X2"  is  the  time  delay  for  transporting  the  fuel  from  the  heat 
exchanger  back  to  the  core. 

Equations  for  the  Core  Pressure,  the  Core  Inlet  and  Exit  Fuel  Mass 
Flow  Rates  and  the  Core  fuel  Mass 

Core  Average  Pressure 

The  core  average  pressure  is  given  by  the  relation 
Pdt)  = n(t)  TKt)  Ru  e / V 

where  "e"  is  a constant  for  conversion  of  pressure  from  atmosphers  to 
pascals  (see  the  end  of  the  Appendix  for  a list  of  symbols  for  the  other 
variables).  This  equation  can  be  written  in  terms  of  the  fuel  mixture 
(UF6  + He)  mass  inside  the  core  as 

Pf  (t)  = ei  Mf(t)  Tf(t)  (B-35) 

where  "ei"  is  a constant. 

For  small  deviations  for  steady  state,  Eq.  (B-35)  can  be  linearized 


as  follows 


5P(t)  = ej  [ff®  5M(t)  + Mf®  8Tf  (t)] 
or,  taking  the  Laplace  transform  of  the  above  equation, 
3P(s)  = C2  3Mf{s)  + C3  9Tf(s) 
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(B-36) 

where  "e2"  and  "ea"  are  constants. 

Core  Fuel  Mass 

The  change  in  the  core  fuel  mass  occurs  due  to  a difference  in 
the  core  inlet  and  exit  fuel  mass  flow  rates,  i.e., 

(dMf/dt)  = m2'(t)  - ms'lt) 

(under  steady  state  conditions,  m2=  013'=  mfss  ) 

Considering  small  variations  from  steady  state,  we  can  write  the  above 
equation  as 

(d3M/ dt)  = 9m2 ' - 9 m3' 

or,  taking  the  Laplace  transforms  of  the  above  equation, 
s 9M({s)  = 9 m2'  (s)  - 9m3'  (s) 
or, 

9Mf  = (1/s)  (9m2’  - 9ni3')  . (B-37) 

(In  Chapter  4.  the  transfer  functions  G13  and  G14  are  given  by 
G13  =l/s  and  G14  = -1/s.) 

Core  Inlet  Fuel  Mass  Flow  Rate 

The  fuel  mass  flow  through  the  inlet  nozzle  is  given  by 
ni2'(t)  = d2'(t)  V2'  (t)  Am 


(B-38) 
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with  d2'  given  by 

d2'(t)  = P(t)  / (RfT2’(t))  . 

In  the  above  equation,  as  an  approximation,  the  core  inlet  pressure  is 
taken  to  be  the  same  as  the  core  average  fuel  temperature,  and 
V2'  (t)  = (yRfT2')i/2  M2'(t)  . 

The  core  inlet  Mach  number,  M2',  can  be  written  in  terms  of  the  core 

average  pressure  and  inlet  stagnation  temperature,  and  after 
simplifying  the  above  eqation,  (B-38),  we  finally  get 
am2'(t)  = fi  3P(t)  + f2  aTo2  (t) 

where 


and  Cg  is  a constant.  In  terms  of  the  Laplace  transform,  we  can  write 
am2'(s)  = fi  ap(s)  + f2  aTo2  (s) . (B-39) 

Core  Exit  Fuel  Mass  Flow  Rate 

At  the  core  exit,  the  fuel  mass  flow  is  assumed  to  always  be 
choked.  So,  we  can  write  the  expression  for  the  maximum  flow 
through  a converging-diverging  nozzle  (which  occurs  when  the  nozzle 
is  choked)  as 
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riig-  (t)  = Pqc  Aex 

VRfToc  • (B-40) 

The  core  stagnation  pressure  is  assumed  to  be  the  same  as  the  core 

average  pressure  and  the  core  stagnation  temperture  is  taken  to  be 

equal  to  the  the  core  average  temperature.  Since  the  fuel  velocity 

within  the  core  is  very  small  (only  around  30  m/s  or  so),  these 

assumptions  should  be  fairly  good. 

Linearizing  the  above  equation,  we  get 

am3'(t)  = gi  aXf+  g2  3P 

where 

gi  = C Pss  (-1/2)  (Tss)-3/2  . 
g2  = C (Tss)(-i/2)  , and 

"C"  is  a constant.  Taking  the  Laplace  transform  of  the  above  equation 
results  in  the  following  equation  for  the  core  exit  fuel  mass  flow  rate 
3m3'(s)  gi  aTf(s)  + g2  aP(s)  . (B-41) 
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List  of  Variables  Used 

The  following  is  a list  of  variables  used  in  this  Appendix: 

Aex  = Throat  area  of  the  converging-diverging  nozzle  at  the  core 
exit  (m2), 

Ahx  = Effective  heat  transfer  area  of  the  heat  exchanger  primary  (m2) 
Arad  = Surface  area  of  the  radiator  (m2). 

Aw  = Wall  area  of  the  core  (m2) , 

CpBe  = Average  specific  heat  at  constant  pressure  of  the  beryllium 

used  as  the  reflector /moderator  (J/Kg-K), 

Cpf  = Average  specific  heat  at  constant  pressure  of  the  fuel  (UFe+He) 

mixture  (J/Kg-K). 

Cphe  = Average  specific  heat  at  constant  pressure  of  the  helium 

used  in  the  secondary  side  of  the  heat  exchanger  (J/Kg-K), 
CpHietai  = Average  specific  heat  at  constant  pressure  of  the  metal  used 

in  the  heat  exchanger  primary/ secondary  side  (J/Kg-K), 
h = Average  convective  heat  transfer  coefficient  for  the  fuel/beiyllium 
wall  heat  exchange  inside  the  core  (w/m2-K), 
hhx  = Average  convective  heat  transfer  coefficient  for  the  heat  transfer 

between  the  primary  and  the  secondary  side  of  the  heat 
exchanger  (w/m2-K), 

M = Mach  number  of  the  fuel, 
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Mbc  = Mass  of  the  beryllium  used  as  the  moderator/reflector  (Kg), 

MfC  = Fuel  mixture  (UFe+He)  mass  within  the  core  (Kg), 

Mf«s  = Steady  state  fuel  mass  inside  the  core  (Kg), 

Mfiix  = Mass  of  the  fuel  within  the  primary  side  of  the  heat 
exchanger  (Kg), 

Mhe  = Mass  of  the  helium  within  the  secondary  side  of  the  heat 
exchanger  (Kg), 

MmetaiP^  = Mass  of  the  metal  used  (e.g.,  as  tubes)  in  the  primary  side 
of  the  heat  exchanger  (Kg), 

Mmetai®®^  = Mass  of  the  metal  used  (e.g.,  as  tubes)  in  the  secondary  side 

of  the  heat  exchanger  (Kg), 
m2-  = Fuel  mass  flow  rate  at  the  core  inlet  (Kg/sec), 

ms-  = Fuel  mass  flow  rate  at  the  core  exit  (Kg/sec), 

mss  = Steady  state  fuel  mass  flow  rate  into  (or  from)  the  core 
(Kg/sec), 

mftix  = Average  fuel  mass  flow  rate  through  the  heat  exchanger 
(Kg/sec), 

N = Core  thermal  power  (watts), 

n = Number  of  moles  of  fuel  mixture  (UFe+He)  within  the  core, 

P = Pressure  variable  (Pascals), 

T2-  = Core  inlet  fuel  temperature  (K), 
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Ta-  = Core  exit  fuel  temperature  (K), 

Tf  = Core  average  fuel  temperature  (K), 

Tbc  = Average  beryllium  temperature  (K), 

Toi  = Fuel  temperature  at  the  heat  exchanger  outlet  (K), 

Tq4  = Fuel  temperature  at  the  heat  exchanger  inlet  (K), 

Thx  = Average  fuel  temperature  in  the  primary  side  of  the 
heat  exchanger  (K), 

Trad  = Average  temperature  of  the  radiator  (K), 
t = time  variable, 

Rf=  Fuel  mixture  gas  constant  (N-m/Kg-K), 

s = Laplace  transform  variable, 

V = Core  volume  (m3) , 

V = velocity  of  the  fuel  (m/sec), 

V = Ratio  of  specific  heats  for  the  fuel  mixture, 

O = Stefan -Boltzmann  constant  (w/m2-K4),  and 
£ = Emissivity  of  the  radiator  material 

Note:  A " - " on  the  top  of  variables  indicate  average  or  mean 


quantities. 
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